{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# 时间序列中最基本的概念\n",
    "很多人看不懂时间序列的公式、很多人看不懂时间序列的数据、代码、数学函数。\n",
    "本文就是将数据、公式、代码、数学函数串联起来，让你不仅仅了解抽象的数学公式，也能将抽象的内容联系到具体的数据和python代码中。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>year</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2001</th>\n",
       "      <td>24123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2002</th>\n",
       "      <td>25159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003</th>\n",
       "      <td>26230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004</th>\n",
       "      <td>27293</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      value\n",
       "year       \n",
       "2001  24123\n",
       "2002  25159\n",
       "2003  26230\n",
       "2004  27293"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"../datasets/近20年城镇就业人员数据(万人).csv\", index_col=['year'])\n",
    "data.head(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(2001, 0, '2001'),\n",
       " Text(2002, 0, '2002'),\n",
       " Text(2003, 0, '2003'),\n",
       " Text(2004, 0, '2004'),\n",
       " Text(2005, 0, '2005'),\n",
       " Text(2006, 0, '2006'),\n",
       " Text(2007, 0, '2007'),\n",
       " Text(2008, 0, '2008'),\n",
       " Text(2009, 0, '2009'),\n",
       " Text(2010, 0, '2010'),\n",
       " Text(2011, 0, '2011'),\n",
       " Text(2012, 0, '2012'),\n",
       " Text(2013, 0, '2013'),\n",
       " Text(2014, 0, '2014'),\n",
       " Text(2015, 0, '2015'),\n",
       " Text(2016, 0, '2016'),\n",
       " Text(2017, 0, '2017'),\n",
       " Text(2018, 0, '2018'),\n",
       " Text(2019, 0, '2019'),\n",
       " Text(2020, 0, '2020')]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAElCAYAAABkuN96AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABKA0lEQVR4nO3dd3hc1bX38e+SZEvuxr1XXLANuAhsSoAQiuk1AUIJaaTekECAEHITQgkkuXlzyU0lIXQwptiY6jRKKAb3hg0xYFuW3I3kJskq6/3jbOFhLMljaWYkjX6f55lHmnPOzN77zJmZNbuauyMiIiIiLUdWU2dARERERA6MAjgRERGRFkYBnIiIiEgLowBOREREpIVRACciIiLSwiiAExEREWlhFMCJSK3M7D4zu62p85FMZnaCma1r6nwkwsx+aGZ/SeLzrTazk5L1fCLStBTAidQhfOGVmtkOMys2szfM7OtmltD7xsyGmJmbWU6q8yqZx91/5u5fachjW2LwbWY3m9lD9ezPNbN7zGxNeE8uMrPT4o75jJmtNLPdZvaSmQ1Ofc5FmoYCOJH6neXunYDBwJ3ADcA9TZul1DOz7KbOg0icHKAAOB7oAvwImG5mQwDMrAfwFPDfQDdgHvBYk+Q00I83SSUFcCIJcPcSd58FXAR8wczGAZjZGWa20My2m1mBmd0c87BXw99iM9tpZkeZ2XAz+5eZbTWzLWb2sJl1rStdMzvazOaaWUn4e3TMvpfN7FYzez3USPwtfInV9jxXmtlrcdvczA4O/99nZn8ws+fNbBfw6XBYDzP7e3j+V2JrNMzsrlDm7WY238w+FbPvZjObbmYPhMcuN7P8eso5OqSzzczeNbPPxey7z8x+b2YvhPP4upn1MbP/NbOPQo3LhJjjV5vZjWb2Tth/r5nl1ZHuIeE8Foc8nh22H2FmG2MDWTM738wWh/+zzOwHZvZ+eC2nm1m3mGOnhBrbYjNbbGYnxL0WH4Tz8qGZXVpH3j6ukYqpzf2Cma0N185NdTzuKuBS4Ppwvp6J2T3ezJaE6+mx2PNiZmeGWq2a2ubDan2x+Pja+U4oxxYz+6WFmun9XeNmdoOZFYbyv2tRrdlU4IfARSHPi+PTdPdd7n6zu69292p3fxb4EJgUDjkfWO7uj7t7GXAzcLiZja4l/9eZ2ZNx235jZneF/7tYVNu3PuT1tpprIYHyrQ5lXALsMrOc2spc17kVSZi766abbrXcgNXASbVsXwt8I/x/AnAo0Y+hw4CNwLlh3xDAgZyYxx4MnAzkAj2Jgrz/rSP9bsBHwOVEtQ+XhPvdw/6XgfeBkUC7cP/OOp7rSuC1uG0OHBz+vw8oAY4JZckL23YAx4X83hX7HMBlQPeQt2uBDUBe2HczUAacDmQDdwBz6shbB6KalS+G55oAbAHGxORtC9EXdR7wL6Iv7ivCc98GvBT3ui0DBoZz+DpwW8zrtS783wZYRRQ4tAVODOUdFfa/A5wW87wzgGvD/1cDc4AB4dz8CXg07OsPbA1lzwqv99bwencAtsek0RcYW8d5uRl4KO5a+nN4rQ8HyoFD6njsfTVljjsvbwP9wnlZAXw97JsAbAImh3P6hXB8bh3P78BL4XkGAe8BX9nfNQ6MCq91v5hyDY8vb4Lvz95E19jocP8u4A9xxywDLqjlsX2BXUDXcD8nlH9SzGv9p/B69Qrn7WuJvIfDeVtEdP21q6/MuunWmJtq4EQOXBHRFxfu/rK7L/WoRmAJ8ChRE0+t3H2Vu//d3cvdfTPw/+o5/gzgP+7+oLtXuvujwErgrJhj7nX399y9FJgOjG9EuZ5299dDWcrCtufc/VV3LwduAo4ys4GhLA+5+9aQt18RfaGNinm+19z9eXevAh4kCjpqcyaw2t3vDc+1EHgS+GzMMTPcfX7I1wygzN0fCM/9GFEAEuu37l7g7tuA24mC33hTgI5EQe8ed/8X8GzMsfcTBamE2rVTgUfCvq8DN7n7unBubgYutKjJ7DLg+VD2anf/O1Fz3unhsdXAODNr5+7r3X15HeelNj9191J3Xwwspu5zWpffuHtROC/PsPd6uQr4k7u/5e5V7n4/UYA4pZ7n+rm7b3P3tcD/Es7bfq7xKqLrZIyZtfGoNu39AywDZtYGeBi4391Xhs0diX6ExCoBOsU/3t3XEwVeNdfYVGCLu883s95Er9V3Par12wT8Grg4gfLV+E24/kqTVWaReArgRA5cf2AbgJlNtqiz9GYzKyH6Yq+1GTMc39vMpoXmlO3AQ/Uc3w9YE7dtTUi/xoaY/3cTfYk1VEF929x9J1G5+wGY2ffNbEVojism6pcUW5b4vOVZ7X2CBgOTQ9NdcXiuS4E+McdsjPm/tJb78eWOLcuamjzH6QcUuHt13LE15/ch4Cwz6wB8Dvh3+OKvyfOMmPyuIPqi7h32fTauPMcCfd19F1Ez/NeB9Wb2XG1NfPVo7Otd1+MHA9fG5XkgtZ+3GrWe4/qucXdfBXyXKODdFI6rL419hKbaB4E9wLdjdu0EOscd3pmoVrU2Hwfo4e+D4f/BRLWz62POxZ+IauISfQ/Hvm8aXWaR2iiAEzkAZnYE0Rd8TX+yR4BZwEB37wL8EbCwz2t5ip+F7Ye6e2eiLw6r5TiIavriR9ENAgobkPVdQPuaO2bWp5ZjasvvwJjHdCSqeSyyqL/b9USBzUHu3pWotqOustSnAHjF3bvG3Dq6+zca8Fz75JvonBXVckwRMNA+Oar44/Pr7oXAm0R9qy5n7xd8TZ5Pi8tzXnhMAfBg3L4O7n5neN7Z7n4yUTPeSqJm0WSr7bWsTwFwe1ye24da37rUdY7rvcbd/RF3P5bo2nbg54nm2cyMaBBRb6Km0YqY3cuJqZEMgffwsL02M4HDLOrPeiZRjR5E56Ic6BFzLjq7+9hEyldbWeops0iDKYATSYCZdTazM4FpRP10loZdnYBt7l5mZkcCn4952Gai5rJhMds6EdUUlJhZf+C6epJ9HhhpZp8PHaEvAsYQNfMdqMXAWDMbHzqu35zg4043s2PNrC1wK1E/toJQjkqiMuaY2Y/Zt/YjUc8SlfNyM2sTbkeY2SENfD6Ab5nZgND0eRO1j0Z8i6gW6vqQ5glEzdPTYo55gChQPZRohGONPwK3WxjUYWY9zeycsK+m5u5UM8s2szyL5p8bEGpvzgnBRTnRtRBbA5gsG/nkdbc/fwa+HmqUzcw6WDRAZ5/mxxjXmdlBoUn9avae4zqvcTMbZWYnmlkuUf+1UvaWfyMwxOqfpucPwCFEo8NL4/bNIGqaviBc4z8GlsQ0sX5CaI5/guhH2NuhKbimefVvwK/C+z7LooELNc2kB/Ie3l+ZRRpMAZxI/Z4xsx1Ev8pvIurv8sWY/d8EbgnH/JioHxoA7r6bqP/V66EpZgrwU2AiUW3Vc3wyKPgEd99KVDNwLVEn+OuBM919y4EWwt3fA24B/gH8h701iPvzCPAToqbTSextcpoNvEjUeX0N0RdTbU2wieRtB3AKUR+jIqJmvp8T9RtqqEeIvoQ/IBrosc+caO6+hyhgO41okMTvgSvivvBnEJpLw+tZ4y6imte/hdd+DtEAAEKAew7R4IjNROflOqLP2yzgmlDObUR9pxpT01iXe4j6XBWb2cz9Hezu84CvAr8lGiizimjgS32eBuYTddh/jr3T69R3jecSTcezheh17gXcGPY9Hv5uNbMF8YmFYPlrRP32Nlg0WnWnhVG8oT/aBUTvuY+IXo+L91OG+4mC8wfjtl9BNLDlnfBcTxDVmO6vfLWpr8wiDWbuB1rTLiLSfJnZaqIRkf9I0vO9TzQCMSnPlwnMzIERoX9Xi2Vmg4iasfu4+/amzo/IgVANnIhIHczsAqI+S/9q6rxIcoWm2muAaQrepCXSLNEiIrUws5eJ+hxeHjdSVVq40AdxI1Hz/9Qmzo5Ig6gJVURERKSFUROqiIiISAvT6ppQe/To4UOGDGnqbIiIiIjs1/z587e4e8/47a0ugBsyZAjz5s1r6myIiIiI7JeZxa/IA6gJVURERKTFUQAnIiIi0sIogBMRERFpYRTAiYiIiLQwCuBEREREWphWNwpVREREpKFmLizkl7Pfpai4lH5d23HdqaM4d0L/tOdDAZyIiIhIAmYuLOTGp5ZSWlEFQGFxKTc+tRQg7UGcmlBFREREEvDL2e9+HLzVKK2o4pez3017XlQDJyIiIlKPqmrnrQ+3UlhcWuv+ojq2p5ICOBEREZE41dXOgrUf8eyS9Ty3dD2bd5RjgNdybL+u7dKdPQVwIiIiIgDuzpJ1JTy7pIjnlqynqKSM3JwsThzdizMP68eu8kp+Mmv5J5pR27XJ5rpTR6U9rwrgREREpNVyd1Zu2MEzi4t4dsl61m7bTZts47gRPbl+6mhOGtObjrl7w6W2OVkahSoiIiLSFFZt2smzS4p4ZnER72/eRXaWcfTw7nz7xIM5dUwfurRvU+vjzp3Qv0kCtngK4ERERKRVWLt1N88siWraVqzfjhlMHtqNLx4zlNPG9aF7x9ymzmLCFMCJiIhIi1fXBLvrS0p5bsl6nllcxOJ1JQBMHNSVn5w1htMP7UvvznlNnPOGMffaxlNkrvz8fJ83b15TZ0NERESSJH6CXYA22caAg9rx4ZbdABzavwtnHtaXMw7ry4CD2jdVVg+Ymc139/z47SmvgTOzbGAeUOjuZ5rZfcDxQEk45Ep3X2RmBtwFnA7sDtsXhOf4AvCjcPxt7n5/2D4JuA9oBzwPXO2tLSIVERFp5WqbYLeiyinYVsr3TxnJmYf1Y0iPDk2Uu9RIRxPq1cAKoHPMtuvc/Ym4404DRoTbZOAPwGQz6wb8BMgnmn5lvpnNcvePwjFfBd4iCuCmAi+ksCwiIiLSTBQWlzJ72YY6J9itqna+feKINOcqPVIawJnZAOAM4Hbgmv0cfg7wQKhBm2NmXc2sL3AC8Hd33xae8+/AVDN7Gejs7nPC9geAc1EAJyIikrFWbdrJ7OUbmL18A0tCn7acLKOyet8GuKaYYDddUl0D97/A9UCnuO23m9mPgX8CP3D3cqA/UBBzzLqwrb7t62rZLiIiIhnC3VletJ0Xl23gxeUbWLVpJwDjB3blB6eN5tSxfVhcULxPH7immmA3XVIWwJnZmcAmd59vZifE7LoR2AC0Be4GbgBuSVU+Ql6uAq4CGDRoUCqTEhERkUaqCstYvbhsAy+GJtLsLGPy0G5ccdRgThnThz5d9o4eHRr6tzWHCXbTJZU1cMcAZ5vZ6UAe0NnMHnL3y8L+cjO7F/h+uF8IDIx5/ICwrZCoGTV2+8th+4Bajt+Hu99NFCySn5+vQQ4iIiLNzJ7Kat78YCuzl2/gb8s3smVnOW2zs/jUiB5cfdIITjqkN906tK3z8c1lgt10SVkA5+43EtW2EWrgvu/ul5lZX3dfH0adngssCw+ZBXzbzKYRDWIoCcfNBn5mZgeF404BbnT3bWa23cymEA1iuAL4v1SVR0RERJKrdE8Vr7y3mdnLN/CPFRvZUVZJh7bZnDC6F1PH9uHTo3t9Yhkr2aspzsrDZtYTMGAR8PWw/XmiKURWEU0j8kWAEKjdCswNx91SM6AB+CZ7pxF5AQ1gEBERaTZqm1z306N78dLKTby4bAMvv7eJsopqurZvw9SxfZg6rg/HHNyDvDbZTZ31Zk8T+YqIiEjS1Ta5bpZFf6sdenfO5dSxfZg6tg9HDu1GTnZWE+W0eWuyiXxFRESk9fn5iyv3mVy32qFjbg4PfPlIxg/oSlZNRCcHTAGciIiIJEVFVTWvvLuZGQsLWV9SVusxu8ormTjooFr3SeIUwImIiEiDuTuLCoqZubCQZ5asZ9uuPXTr0JYOudnsKq/a5/hMnlw3nRTAiYiIyAEr2LabGQsLmbmwkA+27CI3J4uTxvTm/An9OW5kT55bsr7VTa6bTgrgREREJCEluyt4dmkRMxYUMm/NRwBMGdaNrx8/nKmH9qFzXpuPj62Zk601Ta6bTgrgREREpE7llVW8tHIzMxau46WVm9lTVc3BvTp+HIz1r6dJtLVNrptOCuBERETkE9yd+Ws+YsbCQp5dsp6S0gp6dMzlsimDOX9if8b260w0H780FQVwIiIiAsCHW3Z93K9t7bbd5LXJ4tSxfThvQn+OPbiH5mprRhTAiYiItCLxqyN844ThVLvz1IJCFhUUYwZHD+/Odz4zgqnj+mgpq2ZKr4qIiEgrEb86QmFxKT+aGS1JPrpPJ248bTTnjO9Pny55TZlNSYACOBERkVbiZ8+v2Gd1BIBenXJ58bvHNUGOpKEUwImIiGSwktIKZi0uYvrcAjbtKK/1mM11bJfmSwGciIhIhnF35nywjenzCnh+6XrKK6sZ3acTXdrlUFJauc/xWh2h5VEAJyIikiE2lJTx5IJ1TJ9XwJqtu+mUm8OFkwZw0REDObR/F55eVKTVETKEAjgREZEWrKKqmn+u2MT0eQW8/O4mqh0mD+3G1Z8ZwWnj+tKubfbHx2p1hMyhAE5ERKQFWrVpJ9PnFfDUgnVs2bmHXp1y+frxw/lc/kCG9OhQ5+O0OkJmUAAnIiLSQuwqr+S5Jet5bF4B89d8RE6WceLoXlx0xECOH9lTE+22IgrgREREmjF3Z8HaYqbPLeDZJUXs2lPFsJ4duPG00Zw/cQA9O+U2dRalCSiAExERaQbiV0j4+gnDKNtTzfR5Bfxn007atcnmzMP6ctERA5k0+CCtRdrKKYATERFpYrWtkPDfM5cDMGFQV+48/1DOPLyflrWSj+lKEBERaWJ3vFD3CgkzvnlME+RImjsFcCIiIk1gZ3klzy9Zz+PzC9i4XSskyIFRACciIpIm7s7bH27j8fnreH7penbvqWJYjw50zsthe5lWSJDEKYATERFJsaLiUp6cv44nFqxjzdbddGibzdmH9+Oz+QOYOOggrZAgB0wBnIiISAqUVVQxe/kGnpi/jtdWbcEdpgzrxndOHMFph/ahfdu9X8FaIUEOlAI4ERGRJHF3lqwr4fH5BcxaVMT2skr6d23Hf504ggsnDmBQ9/Z1PlYrJMiBUAAnIiLSSJt3lDNzYSGPzy/gvY07yc3JYuq4PnwufyBHDetOVpbmbJPkUgAnIiLSABVV1by0chOPz1/HSys3UVntjB/YldvPG8eZh/WjS7s2TZ1FyWAK4EREROoQvzrCdaeO4pC+nXl8XgEzFxWyZeceenTM5cvHDuXCSQMY0btTU2dZWgkFcCIiIrWobXWE701fhDvkZBmfOaQXn500kONH9aSNFpGXNFMAJyIiUotfzl65z+oI7tClXQ7/uvYEunfUIvLSdBTAiYiIxCjZXcGTC9ZRWFxW6/7tpZUK3qTJKYATEZFWz91ZWFDMI2+t5ZnFRZRXVtMm26io8n2O1eoI0hzsN4Azs97Az4B+7n6amY0BjnL3e1KeOxERkRTaWV7JzIWFPPzWWlas306HttlcMGkAnz9yEKs27dTqCNJsJVIDdx9wL3BTuP8e8BigAE5ERFqkZYUlPPL2Wp5eWMiuPVUc0rczt583jnPG96djbvTVOK5/F0CrI0jzlEgA18Pdp5vZjQDuXmlmVft7kIiISHNSuqeKZ5YU8fBba1lcUExuThZnHd6PSycPYvzArpjtO9muVkeQ5iqRAG6XmXUHHMDMpgAliSZgZtnAPKDQ3c80s6HANKA7MB+43N33mFku8AAwCdgKXOTuq8Nz3Ah8GagCvuPus8P2qcBdQDbwF3e/M9F8iYhI6/CfjTt4+K21PLlgHTvKKjm4V0d+ctYYzp8wgC7tNdmutEyJBHDXALOA4Wb2OtATuPAA0rgaWAF0Dvd/Dvza3aeZ2R+JArM/hL8fufvBZnZxOO6i0OfuYmAs0A/4h5mNDM/1O+BkYB0w18xmufs7B5A3ERHJQOWVVby4bAMPz1nL26u30TY7Wtrq0smDOHJot1pr20Rakv0GcO6+wMyOB0YBBrzr7hWJPLmZDQDOAG4HrrHoHXMi8PlwyP3AzUQB3Dnhf4AngN+G488Bprl7OfChma0CjgzHrXL3D0Ja08KxCuBERFqp1Vt28ejba3l8/jq27drD4O7tufG00Vw4aYCm/pCMUmcAZ2bn17FrpJnh7k8l8Pz/C1wP1Kwt0h0odvfKcH8dUNO5oD9QAB/3sysJx/cH5sQ8Z+xjCuK2T66jLFcBVwEMGjQogWyLiEhzFb+81TUnj6B92xwefmstr63aQnaWcfIhvbl0yiCOGd5DC8lLRqqvBu6sevY5UG8AZ2ZnApvcfb6ZnXDgWUsed78buBsgPz9/30l9RESkRahteatrH18CQL8ueVx78kg+d8RAenfOa8psiqRcnQGcu3+xkc99DHC2mZ0O5BH1gbsL6GpmOaEWbgBQGI4vBAYC68wsB+hCNJihZnuN2MfUtV1ERDLQL2pZ3gqge4e2/PuGE8lWbZu0EvtdfdfMupvZb8xsgZnNN7O7wqjUern7je4+wN2HEA1C+Je7Xwq8xN5BEF8Ang7/zwr3Cfv/5e4etl9sZrlhBOsI4G1gLjDCzIaaWduQxqwEyy0iIi3I6i27+PmLKymqY3mrbbv2KHiTViWRUajTgFeBC8L9S4km8j2pgWneAEwzs9uAheydEPge4MEwSGEbUUCGuy83s+lEgxMqgW+5exWAmX0bmE00jchf3X15A/MkIiLNTFlFFbOXb2Da2wW8+cFWsgzycrIoq6ze51gtbyWtjUWVXPUcYLbM3cfFbVvq7oemNGcpkp+f7/PmzWvqbIiISB3e27iDR99ey4yFhRTvrmBgt3ZclD+QCycNZM4HW2td3uqO8w/VhLuSkcxsvrvnx29PpAbub2Fetunh/oVEtV4iIiJJsau8kueWrOfRuWtZuLaYNtnGKWP7cMkRgzh6ePePR5LWBGla3kpau0Rq4HYAHYCaOussYFf43929c60PbKZUAyci0jy4O0sLS3j07QKeWVzEzvJKhvfswCVHDuK8Cf01b5sIjaiBc/dO+ztGREQkUSWlFTy9qJBH3y5gxfrt5LXJ4oxD+3HxkQPJH3yQVkkQSUAiTaiY2dnAceHuy+7+bOqyJCIimcbdmbv6I6a9vZbnlq6nvLKaMX07c+s5Yzl7fH+6tNOapCIHYr8BnJndCRwBPBw2XW1mx7j7jSnNmYiItBjxqyPU9EvburOcJxesY9rcAj7YvIuOuTlcMGkAlxwxiEMHdGnqbIu0WIn0gVsCjHf36nA/G1jo7oelIX9Jpz5wIiLJFb86AkDb7CwO6duRd9bvoKLKmTioKxcfOYgzD+tL+7YJNf6ICI0bhQrQlWhuNohWSBAREQGiEaHxqyPsqapmSeF2vnj0UC4+ciAje6s7tUgyJRLA3QEsNLOXACPqC/eDlOZKRERaBHensLi0jp3w47PGpDdDIq1EIqNQHzWzl4n6wQHc4O4bUporERFp1naUVfDUgkIenLOmzmO0OoJI6iQyiMGAzwDD3P0WMxtkZke6+9upz56IiDQnK9Zv56E5a5ixsJDde6o4tH8XLj5yIDMXFlJWsXeJq3Ztsrnu1FFNmFORzJZIE+rviSbxPRG4BdgBPMneGjkREclgeyqreWHZeh6as4a5qz8iNyeLsw7vx+VTBnP4wK4ATBnaXasjiKRRIgHcZHefaGYLAdz9IzNrm+J8iYhIEyssLuWRt9bw2NwCtuzcw+Du7bnp9EO4cNIADurwya+Bcyf0V8AmkkaJBHAVYeoQBzCznuxdVktERDJIdbXz2qotPDhnDf9csREHPjO6F5dNGcxxI3p+vCapiDStRAK43wAzgN5mdjvRYvY/SmmuREQkrYp37+GJ+et4aM4aVm/dTfcObfn68cO55MhBDOzWvqmzJyJxEhmF+rCZzScayABwrruvSG22REQkHZauK+GBN1cza3ER5ZXV5A8+iO+dPJKp4/qQm5Pd1NkTkTokOpFve6CmGVXjwkVEWojalriaOq4Pzy5Zz4Nz1rC4oJj2bbO5YNIALps8mDH9Ojd1lkUkAYkspfVj4LNEI08NOBd43N1vS3nuUkBLaYlIa1HbElc5WUabbKO0oprhPTtw+ZTBnD9pAJ3ztJi8SHPUmKW0LgUOd/ey8ER3AouAFhnAiYi0FrUtcVVZ7eRkG498dTJHDetONNWniLQ0WQkcUwTkxdzPBQpTkx0REUmGktKKOpe4Kq+o5ujhPRS8ibRgidTAlQDLzezvRH3gTgbeNrPfALj7d1KYPxEROQCrNu3g/jfW8OSCdXUeoyWuRFq+RAK4GeFW4+XUZEVERBqiutp56d1N3PfGav79ny20zc7i7PH9GNy9Pb9/6f1PNKNqiSuRzJDINCL3pyMjIiJyYLaXVfD4vHU88OZq1mzdTe/OuXz/lJFcfOQgenTMBWDgQe21xJVIBkp0GhEREWkmVm3ayQNvruaJ+evYvaeK/MEH8f1ToulB2mR/smuzlrgSyUwK4EREWoDqaueV9zZz7xurefW9zbTNjhaUv/LoIRw6oEtTZ09E0kwBnIhIM7a9rIInQjPp6q276dUpl2tPHsklk/c2k4pI61NnAGdmzxAWsK+Nu5+dkhyJiAjvb97JA29EzaS79lQxafBBXHPKKE6rpZlURFqf+mrg/if8PR/oAzwU7l8CbExlpkREWqOaZtL73ljNK6GZ9MzD+3Ll0UM4bEDXps6eiDQjdQZw7v4KgJn9Km4Jh2fMTGtRiYg0UPz6pP914nBKK6p54M01fLhlF7065XLNySO55MhB9OykZlIR2VcifeA6mNkwd/8AwMyGAh1Smy0RkcwUvz5pYXEpP3hqGQATB3Xle5dMYOrYPrTNUTOpiNQtkQDue8DLZvYB0WL2g4GvpTRXIiIZ6pezV+6zPilAz465PPXNY5ogRyLSEiUyke+LZjYCGB02rXT38tRmS0Qks5RVVDFzYSGFxWW17t+yUx+rIpK4RKcRmQQMCccfbma4+wMpy5WISIbYtL2MB+es4eG31rJt1x5ysozK6n0H+Gt9UhE5EPsN4MzsQWA4sAioqfd3QAGciEgdlhWW8NfXP+SZxUVUVjufGd2bLx87lA0lpfxwxjKtTyoijZJIDVw+MMbd65wTTkREoKra+ceKjfz1tQ9568NttG+bzaWTB3Pl0UMY0mPv2C8z0/qkItIoiQRwy4jmgVuf4ryIiLRIO8sreXxeAfe9ES0q379rO246/RA+d8RAurRrs8/xWp9URBorkQCuB/COmb0NfNzLVisxiEhrV7BtN/e/sZrH5hawo7ySSYMP4oapozllTG9ytFqCiKRQIgHczQ15YjPLA14FckM6T7j7T8zsPuB4oCQceqW7LzIzA+4CTgd2h+0LwnN9AfhROP42d78/bJ8E3Ae0A54HrlZTr4ikkrszf81H3PPah8xevoEsM04/tC9fOnYo4wd2bersiUgrkcg0Iq808LnLgRPdfaeZtQFeM7MXwr7r3P2JuONPA0aE22TgD8BkM+sG/ISoL54D881slrt/FI75KvAWUQA3FXgBEZEkq6iq5vml67nntQ9Zsq6ELu3a8LXjh3PFUYPp20UjSEUkvRIZhToF+D/gEKAtkA3scvfO9T0u1ITtDHfbhFt9tWPnAA+Ex80xs65m1hc4Afi7u28L+fk7MNXMXgY6u/ucsP0B4FwUwIlIEhXv3sPDb63lgTdXs3F7OcN6duC2c8dx/sT+tG+b6ExMIiLJlcinz2+Bi4HHiWrBrgBGJvLkZpYNzAcOBn7n7m+Z2TeA283sx8A/gR+EiYH7AwUxD18XttW3fV0t22vLx1XAVQCDBg1KJOsi0srEr096xVGDWbttN08uWEdZRTWfGtGDO88/jONH9iQry5o6uyLSyiX089HdV5lZtrtXAfea2ULgxgQeVwWMN7OuwAwzGxcet4GoNu9u4AbglgbmPyHufndIi/z8fPWRE5FPqG190jteWEm2wYWTBvKlY4cyqk+nJs6liMheiQyT2m1mbYFFZvYLM/tego/7mLsXAy8BU919vUfKgXuBI8NhhcDAmIcNCNvq2z6glu0iIgfkFy/WsT5p5zx+fuFhCt5EpNlJJBC7PBz3bWAXUTB1wf4eZGY9Q80bZtYOOBlYGfq1EUadnks0zxzALOAKi0wBStx9PTAbOMXMDjKzg4BTgNlh33YzmxKe6wrg6cSKLSIC28squPvV9ykqqX190o11bBcRaWqJjEJdE/4tM7PfAAPdfVUCz90XuD/0g8sCprv7s2b2LzPrCRjR8lxfD8c/TzSFyCqiaUS+GNLfZma3AnPDcbfUDGgAvsneaUReQAMYRCQB60tKuff11Tz61lp2lFeSm5NFeWX1PsdpfVIRaa4SGYX6MnB2OHY+sMnMXnf3a+p7nLsvASbUsv3EOo534Ft17Psr8Ndats8Dxu2nCCIiAKzcsJ27X/2AWYuKcOD0Q/ty1aeG8f7mnZ/oAwdan1REmrdEBjF0cfftZvYVomk+fmJmS1KdMRGRZHB33nx/K3969QNeeW8z7dtmc9mUwXz52KEM7NYegEMHdAHQ+qQi0mIkEsDlhH5rnwNuSnF+RESSorKqmueWrufP//6AZYXb6dExl+tOHcWlkwfRtX3bfY7X+qQi0pIkEsDdQjSQ4DV3n2tmw4D/pDZbIiINs6u8ksfmFnDPax9SWFzKsJ4duPP8Qzl3Qn/y2mQ3dfZERJIikUEMjxNN4ltz/wMSGIUqIpJOm3aUcf8bq3lozlpKSis4ckg3fnr2WE4c3UsT74pIxklkEMO91LIElrt/KSU5EhE5AKs27eTPr37AjIWFVFRXM3VsH646bhgTBh3U1FkTEUmZRJpQn435Pw84DyhKTXZERPbP3Zm7+iPufvV9/rFiE7k5WXzuiAF85dhhDOnRoamzJyKScok0oT4Ze9/MHgVeS1mORERixK5R2rdrHqeM6c2ighIWFRTTrUNbvnvSCC6fMpjuHXObOqsiImmT0FqocUYAvZKdERGRePFrlBYVl3HfG2vo3qENt547jgsnDqBdWw1MEJHWJ5E+cDuI+sBZ+LuBaAF6EZGUuuOFFbWuUZrbJpvLpwxughyJiDQPiTShahVnEUmbiqpq/rliIw+/tZaN28trPWZ9sdYoFZHWLaEmVDM7Gzgu3H3Z3Z+t73gRkQNVVFzKtLkFPDY3Ctz6dcmjU14OO8oq9zlWa5SKSGuXSBPqncARwMNh09VmdrS7/zClORORjFdd7bz6n808/NZa/rliIw6cMLInt587mE+P7sUzi4u0RqmISC0SqYE7HRjv7tUAZnY/sBBQACciDbJlZzmPz1vHI2+voWBbKT06tuXrxw/nkiMHfbw+KfDx0lZao1RE5JMSHYXaFdgW/u+SmqyISCZzd976cBsPv7WWF5etp6LKmTKsG9efOppTx/ahbU5WrY/TGqUiIvtKJIC7A1hoZi8RjUQ9DvhBSnMlIhmjpLSCpxas4+G31rJq00465+Vw2ZTBXDp5MAf36tjU2RMRaZESGYX6qJm9TNQPDuAGd9+Q0lyJSIvm7ixZV8LDb61h1uIiyiqqOXxgV3554WGceVg/zd0mItJIiTahZgFbwvEjzWyku7+aumyJSEu0e08lTy8q4uG31rCscDvt22Zz3oQBXDp5EOP6q/eFiEiyJDIK9efARcByoDpsdkABnEgrFbu8Vb+u7bhsyiDWl5QxY0EhO8orGdW7E7eeM5ZzJvSnc16bps6uiEjGSaQG7lxglLvXPqOmiLQq8ctbFRaX8vMX3yXb4Ozx/bl08iAmDT4IM2vinIqIZK5EArgPgDaAAjgR4WfP1768Va/Oefz6ovHpz5CISCtUZwBnZv9H1FS6G1hkZv8kJohz9++kPnsi0hyU7K5g1uJCps9bx6Ydtf+W21Ci5a1ERNKlvhq4eeHvfGBW3D5PTXZEpLmornbe/GAr0+cV8OKyDZRXVnNI3850aZdDSamWtxIRaUp1BnDufj+AmV3t7nfF7jOzq1OdMRFpGoXFpTwxbx2Pzy9g3UeldM7L4aIjBvK5/IGM7deZpxdpeSsRkaaWSB+4LwB3xW27spZtItJClVdW8bflG5k+r4DXVm3BHY45uDvXnTqKU8f2Ia/N3nnbtLyViEjTq68P3CXA54GhZhbbhNqJvctqiUgLtryohMfnrWPmokKKd1fQv2s7vnPiCC6cNOATa5LG0/JWIiJNq74auDeA9UAP4Fcx23cAS1KZKRFJnZLdFTy9uJDp8wpYVridttlZnDK2NxcdMZCjh/cgO0vTf4iINHf19YFbA6wBjkpfdkQkFaqrnTfeDwMSlm9gT2U1Y/p25qdnj+Wc8f3o2r5tU2dRREQOQKJLaYlIMxe/OsJ1p44if8hBPDF/HY/PW0dhcTQg4eIwIEFLW4mItFwK4EQyQG2rI1wzfRHVYcKfYw/uwQ2njeaUMb0/MSBBRERaJgVwIhngl7Pf3Wd1hGqHTnk5PP+dT9U7IEFERFqeRBazPwa4GRgcjjfA3X1YarMmIolYsPYjCotLa923s6xSwZuISAZKpAbuHuB7RCsy7LsAooikXWVVNS8u38A9r33IwrXF0a+qWo7T6ggiIpkpkQCuxN1fSHlORGS/tpdV8NjbBdz3xmoKi0sZ3L09Pz17LHk5Wdz8zDtaHUFEpJVIJIB7ycx+CTzFJxezX5CyXInIJ6zZuot7X1/N4/MK2LWnislDu/GTs8bwmUN6fzxvW26bbK2OICLSSiQSwE0Of/NjtjlwYvKzIyI13J25qz/iL//+gL+v2Ei2GWcd3o8vHzu01ilAtDqCiEjrsd8Azt0/3ZAnNrM84FUgN6TzhLv/xMyGAtOA7kT96i539z1mlgs8AEwCtgIXufvq8Fw3Al8m6oP3HXefHbZPJVqTNRv4i7vf2ZC8ijQnFVXVPLdkPfe89iFLC0vo2r4N3zxhOFccNYTenfOaOnsiItIMJDIKtTfwM6Cfu59mZmOAo9z9nv08tBw40d13mlkb4DUzewG4Bvi1u08zsz8SBWZ/CH8/cveDzexi4OfARSG9i4GxQD/gH2Y2MqTxO+BkYB0w18xmufs7B3YKRJqH4t17eOTttTzwxho2bC9jWM8O3H7eOM6fMIB2bTV3m4iI7JVIE+p9wL3ATeH+e8BjRKNT6+TuDuwMd9uEW03T6+fD9vuJpij5A3BO+B/gCeC3ZmZh+zR3Lwc+NLNVwJHhuFXu/gGAmU0LxyqAkxblg807+evrH/Lk/EJKK6o49uAe3HH+oRw/sidZWpdURERqUWcAZ2Y57l4J9HD36aEZE3evNLOEphMxs2yiZtKDiWrL3geKw/NCVHNW02mnP1AQk0YJUTNrf2BOzNPGPqYgbvtkamFmVwFXAQwaNCiRrIuklLvz5vtb+ctrH/KvlZtom53FuRP68aVjhzK6T+emzp6IiDRz9dXAvQ1MBHaZWXfCNFNmNgUoSeTJ3b0KGG9mXYEZwOhG5baB3P1u4G6A/Pz82qbLEkmZ2DVK+3bJ4/hRPVm4tpiVG3bQvUNbrv7MCC6bMpienXKbOqsiItJC1BfA1bTdXAPMAoab2etAT+DCA0nE3YvN7CXgKKBrTO3eAKAwHFYIDATWmVkO0IVoMEPN9hqxj6lru0izEL9GaVFJGY++XUCfzrn8/IJDOWd8f61NKiIiByyrnn09zewa4ASi2rNfAC8AfwZO2t8Tm1nPUPOGmbUjGmywAniJvQHgF4Cnw/+zwn3C/n+FfnSzgIvNLDeMYB1BVDs4FxhhZkPNrC3RQIdZCZRZJG1uf27FPmuUAmRnGRcdMUjBm4iINEh9NXDZQEf21sTVSHRhxb7A/aEfXBYw3d2fNbN3gGlmdhuwkL2DIe4BHgyDFLYRBWS4+3Izm040OKES+FZomsXMvg3MDnn9q7svTzBvIilTuqeKZ5YU8dCcNWzeWV7rMUXFZWnOlYiIZBKLKrlq2WG2wN0npjk/KZefn+/z5s1r6mxIBvpwyy4enrOGx+evo6S0ghG9OrJpRzklpRX7HNu/azte/4HmwhYRkfqZ2Xx3z4/fnkgfOBGpQ1W186+Vm3hwzhpefW8zOVnGqeP6cPmUwUwe2o2nFxV9og8caI1SERFpvPoCuM+kLRciLczmHeVMn1fAI2+tpbC4lD6d87jm5JFcfMRAesWsllCztJXWKBURkWSqM4Bz923pzIhIc+fuzF/zEQ/OWcPzS9dTUeUcc3B3/vvMQzjpkN7kZNc+JkhrlIqISLIlshKDSKu2q7ySmYsKefDNNazcsINOeTlcNmUwl04ezMG9OjZ19kREpBVSACdSh1WbdvDgm2t4ckEhO8srGdO3M3eefyhnj+9H+7Z664iISNPRt5BIjIqqav7+zkYefHMNb36wlbbZWZxxWF8umzKYiYO6Ei3PKyIi0rQUwEmrFLu8Vb+u7bjquKFs21XBo2+vZdOOcvp3bcf1U0dxUf5AunfUElciItK8KICTVid+eavC4lJ+MusdAE4Y1ZM7pgzmhFG9yM5SbZuIiDRPCuCk1bnjhdqXt+rdOZf7vnhkE+RIRETkwCiAk1ZhR1kFLyzbwIwFhWzcXvvyVpvq2C4iItLcKICTjFVV7by2agtPLVjH7OUbKKuoZkj39nTKy2FHWeU+x/fr2q4JcikiInLgFMBJxlmxfjszFhYyc2Ehm3aU06VdGy6YOIDzJw5g4qCuWt5KRERaPAVwkhE27Shj1qIinlxQyIr128nJMj49uhcXTOzPp0f3Ijcn++NjtbyViIi0dArgpMUq3VPF397ZwIyFhbz63maqHQ4f2JVbzhnLmYf1o1uHtnU+VstbiYhIS6YATlqU6mrn7dXbeGrBOp5fuoGd5ZX065LHN04YznkTBmhpKxERaRUUwEmL8MHmncxYWMhTCwopLC6lQ9tsTj+0L+dPHMDkod3I0pxtIiLSiiiAk2YjfnWEb54wnGp3nlxQyKKCYrIMjh3Rk+unjuKUMX1o1zZ7/08qIiKSgRTASbNQ2+oIN81cBsDoPp246fRDOGd8P3p1zmvKbIqIiDQLCuCkyVVXO7c9906tqyP06pTLi989rglyJSIi0nwpgJMms3LDdmYuLGLWokK27NxT6zGbd2h1BBERkXgK4CStCotLeXpRIbMWFbFyww6ys4zjRvSgtKKKj3ZX7HO8VkcQERHZlwI4Sbni3Xt4bul6nl5YxNurtwEwcVA0X9sZh/ale8fcffrAgVZHEBERqYsCOEmJsooq/rFiIzMXFvHKe5uoqHKG9+zAtSeP5Jzx/RnUvf0njtfqCCIiIolTACdJU1XtvPH+FmYuLGL28miS3V6dcvnCUUM4d0J/xvbrjFnd87VpdQQREZHEKICTRnF3lhaWMHNhEc8sKWLzjnI65eZw2rg+nDehP5OHdSdbk+yKiIgklQI4qVf85Lo1zZqrt+zi6UVFPL2okA+27KJtdhafHt2Tc8dHi8fntdEkuyIiIqli7t7UeUir/Px8nzdvXlNno0WobWBBm2yjb5c81m4rxQwmD+3GueP7c9q4vnRp36YJcysiIpJ5zGy+u+fHb1cNnNTpl7Pf3Wdy3Yoqp6i4jBtPG81Zh/fTNB8iIiJNQAGc7KOyqpp/r9pCYXFprfurqp2vHT88zbkSERGRGgrgBIgGIyxeV8LMhYU8u6SILTv3YAa1tbCr1k1ERKRpKYBr5VZv2cXMRYU8vaiID7fsom1OFicd0otzxvdnR2kF//30ck2uKyIi0swogGuFtu4s59kl65m5qJCFa4sxgylDu/ON44dz6rg+dGm3dzBCTnaWJtcVERFpZhTAtRKle6r42zsbeHpREa+8t5mqamd0n07ceNpozh7fj75dam8W1eS6IiIizY8CuAxWWVXNG+9vZeaiQmYv28CuPVX07ZLHVz81jHMn9GN0n85NnUURERFpAAVwGcbdWVa4nZmLCpm1OKyMkJfDWYf345zx/Zk8tBtZWhlBRESkRVMA10LFr5DwpWOHULqnihkLC3l/896VEc6b0J8TRmllBBERkUySsgDOzAYCDwC9AQfudve7zOxm4KvA5nDoD939+fCYG4EvA1XAd9x9dtg+FbgLyAb+4u53hu1DgWlAd2A+cLm770lVmZqL+BUSCotLufXZFQAcObQbX/nUME7XyggiIiIZK5U1cJXAte6+wMw6AfPN7O9h36/d/X9iDzazMcDFwFigH/APMxsZdv8OOBlYB8w1s1nu/g7w8/Bc08zsj0TB3x9SWKYmt7O8kp8+s3yfFRIAenfOZfrXjmqCXImIiEg6pSyAc/f1wPrw/w4zWwHUN5zxHGCau5cDH5rZKuDIsG+Vu38AYGbTgHPC850IfD4ccz9wMxkYwJVVVPHyu5t4ZvF6/rlyI2UV1bUet2l7eZpzJiIiIk0hLX3gzGwIMAF4CzgG+LaZXQHMI6ql+4gouJsT87B17A34CuK2TyZqNi1298pajo9P/yrgKoBBgwYloUSpV1FVzeurtjBrcRF/W76RneWVdO/Qls/lD+T5pevZsnPflmKtkCAiItI6pDyAM7OOwJPAd919u5n9AbiVqF/crcCvgC+lMg/ufjdwN0B+fn4ti0M1D9XVztzV23hmSRHPL93Atl176JSXw2nj+nD2+H4cNaw7OdlZTBx00Cf6wIFWSBAREWlNUhrAmVkbouDtYXd/CsDdN8bs/zPwbLhbCAyMefiAsI06tm8FuppZTqiFiz2+xXB3lhaWMGtREc8uWc+G7WXktcnipEN6c/bh/Th+VE9ycz45grRmYl2tkCAiItI6pXIUqgH3ACvc/f/FbO8b+scBnAcsC//PAh4xs/9HNIhhBPA2YMCIMOK0kGigw+fd3c3sJeBCopGoXwCeTlV5ku0/G3cwa3ERzywuYvXW3bTJNo4f2ZMbTx/NSYf0pkNu/S+NVkgQERFpvVJZA3cMcDmw1MwWhW0/BC4xs/FETairga8BuPtyM5sOvEM0gvVb7l4FYGbfBmYTTSPyV3dfHp7vBmCamd0GLCQKGJutgm27Pw7aVm7YQZbBUcO7840ThnPq2D50bd+2qbMoIiIiLYC5N9suYSmRn5/v8+bNS8lzx0+ue92pozh6eHeeXbKeZ5YUsXBtMQATB3Xl7MP7cfphfenVKS8leREREZGWz8zmu3t+/HatxJAktU2ue830RVSH+PiQvp25YepozjysLwO7tW/CnIqIiEhLpwAuSX45+919JtetduiUl8OMbx7Nwb06NVHOREREJNNkNXUGMkVRcWmt23eWVSp4ExERkaRSAJckdU2iq8l1RUREJNkUwCXJdaeOol2bT87Xpsl1RUREJBXUBy5JNLmuiIiIpIsCuCTS5LoiIiKSDmpCFREREWlhFMCJiIiItDAK4ERERERaGAVwIiIiIi2MAjgRERGRFqbVLWZvZpuBNSlOpgewJcVppDutTCxTpqaViWVKZ1qZWKZ0ppWJZcrUtDKxTOlMK13pDHb3nvEbW10Alw5mNs/d8zMprUwsU6amlYllSmdamVimdKaViWXK1LQysUzpTCudZaqNmlBFREREWhgFcCIiIiItjAK41Lg7A9PKxDJlalqZWKZ0ppWJZUpnWplYpkxNKxPLlM600lmmfagPnIiIiEgLoxo4ERERkRZGAZyIiIhIC6MAThrMzCxD08qo94WZ5aY5vbS9VumSiWWC9JUrU9+/6UgrzeXpmK60QnoZ9x2SzjJl1BdVc5amN3qembVNQzqdzaydp6EDpZkNM7PO7u6pfmOY2ZFm1sPdq1OZTkjrRDObkoYyfRr4avg/pdegmQ0ys4OAnFSmE9Jql6ZrvbuZdUjHtR6TZnYa0uhqZu1TXS4z62NmnVKZRkxaY82su7tXp+FaP87M+qT6s8LMTgLONrO8VKYT0vo08N9mlpOG8zfKzPoTTYSbUmbWLU3Xeto/KxTApVB4k38JINUfKmZ2JvBX4EUzOyFsS3p6IZ2HgFlmdnKynz8urYnA68CNZtYtlUGcmZ0C3AP0T8Xzx6U1FXgA6FLzZk9FuUKZngR+ZWYDUvllY2ZnA9OAx4BLw7ZUvVbnEL1Wj5rZKWY2OEXpnA88CjxnZl81s8mpSCekdYqZ3Qjg7lUp/qw4i+g9/IKZfT5VAZaZnQE8QnRdfNHMslN4TYwBXgJ+a2a9U/l5G95X9wMpue5i0jkVuA/Y5e5lYVuqzt+pRK/VN4GDw/lLVVpnAQ8DvwW+ZmbtU/hanUP0Wj1kZpebWUom3U3nZ8UnuLtuKbgBpwAfAa8A18Zsz0pBWlOBpcCngG8B84FRKUjnTGABMIWoVudNoEMKz2FXYDZwB3An0DNF6ZwdyjUx3M9OUToGdAT+DpwWtnUAcoG8FL1WY4Abw/lrm6JyjQeWAYeGdP8GdEpRWoeHa/0w4DyiL53/AcYkOZ1+wLvAxPBe/iHwR+DkFJTpOGATsBL4n5jtqfisODm8VvnAZ4HngckpSOcMYCFwBHA68C/goFRcEyG9HKJA8S5gOjAgRemcCiwGpoT7ucl+ncLnRB7Rj68LwrYu4Zb0z0DgrPBZMRz4L+AJoGOKzt+omM+KicAMoHuKrvWR4bNiTHiP/ZIoyDo2yemk7bMi/qYauNQ5GPgF8F3gKDO7FpJfExeq1k8FbnX3f7v774AXgdOSlUZIpw1wAvBDd58DvADsAW42s5PNrF+S08si+lDeA3xI9KH2FTM72swmJDMt4EKgv7svMLPOwC/M7AEzOy+Z5fLo3V4GrAdmm1lXoi+dh4FbQnNJo4VmzPOB69z9HaLAYDCQHfYn+5f1IOAdd18KvAp0Bn5jZt80s/FJTmtwSGuJu88AXgYmA2eZ2T5rBTZCNrDW3Re4+9+IXqfFwPlmNimJ6UD0BXATcAwwwcx+BR9/ViStOTU81zFEQeI8d3+c6PX6bNifzOtiEvDf7j6XKDjoAvzczC41s8OSmE7NZ0VNE+PLwHKi99MpZnZ8MtMCTgLaufuccL39lqgm+FvJKpdHyojW7J5jUb+0mURzjv2vmV2SjHTg4++Pk4Dr3f194DVgJ9A77E92jNADWBc+K94jChp/B9yWrM+/WtJ6x91fJSpbN+CiUGObLDmk77PiExTApYi7/x74PdGvjd8CU8zsurCvOll9GsIb/bdETac1H/bbgRE1xyTjg9ndK4AfuPuLocnlSWAO0YV6GdGHQNK4e7W7byEKFN8mepPX/GLrCcn7wnH3K4DXzGwl8AywFpgLXAB8OhlpxKRVCVQRPoyBx4lqx9YQBSFdklCuYuDb7v7PkOYMoA/RL9CaQDKZ3gZ6mNl0YAUwC3gKGAicluSms6VAhZldEe73CWkeDgxLUhq4ewGwzcz+J9z/gKhmcSNR7UEyr79pwBPuvhX4MnC4mf067KsKgX4y0qki+kyaYQFQBPQK+92SNODF3W9x92fNrD3Re/Y5oiBkHHCmmWUl8fxVu/tOos+Kcnf/KVFA9zihj1WyAhF3vw54xczmAk8T1TI+RfQjZmrMeW2wuMf/BbgVuBe4hqhM5yer20D4/rje3f8RNi0iahW4JexPdreLuUA7M/sH8D5RU/7/I7oOz7Wob2ayPiuWAdvN7Efh/gSimrIyotcrKdx9bUgn5Z8VtSWuW/KqUk8APgdcFrc9jygQeBz4ItEv3stpRLVxSOsi4HO17DsD+EX4/yLglCSk8/m47aNj/r+YKPBpk6Tz9/mYbTcSNdceBxQQfRncBnRPUlpfiNn2KHBHzP1LiIKRnCSldXm4Pxj4NdGHWW7YNjKUrcHNTDGv1SUx27LC3zFEzY2HNPY6jyvTJTFlOpsoEKk55jiiJrr2SUrr4nD/K0TNPC8Az4Vt3wTuamQ6A4j6JdbcH08U8Hw/ZtsZRDXcjWryjk8rbt8w4B/Aj4lqh29ozHsrNq34zxyi2st7wv+XAFfQwC4EcelY7PaY/z8dPityk3X+Yq7xr4ZzdQxRrf1j4Trpm+zXCvg/4KaY+ycCzzamXLVcf9lEP1zfA/qFbV2IPqcGJqFMXeO21ZzH7uF92+Dvjf2UKw84FngwZttoos/abkm8LnKIugzMDM/9dNh+CVG/OGtEOscQtXJcFu6PAf6Uis+K+m6qgUsSi0bwPEoU2X/XzH5v0SgbPPqV8wbRB/INRL+qFnoDf93EpDUQuMHM/hDX1LcHqDKzzxP9klqThHSuqUnHzMzdV8b9qiihETW6cefvGjP7Y2gKnEnUt+9h4DvATwAnNAcmIa3/MrM/m1lXd7/E3W+MKVcWUW1mssr1XTP7bXi+F4Bqoi8CiD4AGjyCM+61ujZcf/1irrGt4bmPbWhZaklrEHCdmf0OKHX3WcCGcN1B9GWTAzR4tGhcWteHX7n/AL5EdD2cGw7NIboGG5rOueF5vxzTFFtTIzvczO4K2zoCFTTu+otNq6aG6OP3kke/4E8Hvgb8GXjeoxrwRqdVy2dOFVBtZlcSfT697VFNXWPT8ZhatsKYQ7sTvX/bHHBh6k6rpkwziWpZniKqrbqSqIUgqa8VgLv/F1HNeY0eROeyQeWqI51q4DdE/SPvCefyVGAIUNmQdOLS+lJcmWoGLuwmqokb39A0aknr4/eVu5e5+2vATjP7QTj0YKI+z8l6rXq6e6W7/50o0Loq/AXoBBR7iLIakM7pRC0o44AfmtnPPOqqMoMkf1bsV6oiw9Z0I+qf9Qvge+F+HlHA8RugT8xxXyNqnmtwh+tE0iJ6k+8iChobVOOyn3R6xxz3FWAeMC4FZfo1MJaoY/JpMcc2eOBEPWn9Nq5cXw3lOjTJadV0uh9AVOv2IFFNwTzgsBRffxcR9Q9qRwN/fe7n/PUP18NTRAHqkoaWqZ60Hg3l6hFz3PeJvqjHNjCdnsA/iZqpbiXqyN0zJs2x4TV6nqgJd0IjyhSf1rdjyxJz3IXA6oaWaX9phXNrRE3PO4g+K0anuEzfJBpg1Zj3VH1l6kj0A/n4mDI2eODOAZTrWzTis6KOdHrG7M8jChj+RDQqPyXnL+64qcAHRM2pDf2sqDctog7/s4hGDzf2s6LO93Dccd8NnxUNfa1GhNf62HB/CNGPvM5En6tJ+6xIKD+pfPLWdCNqRvwjIQgA2hN1Zvx93MXT4DdfomkRdUD9F4348E8wnVFEtYmpKtN0PtmkmZTRoQmUq1/Yn6pyPQb8X8wx3YHOqb7+wvZGjwSsp0y/DPeHEo1sG5Tq1ypsuwU4vBFptCXqq5JL1O/xLqIvgD5xx/WmjmbPJKTVK+yvacb6Ao1s7t5fWjHn87nGXOsJlCknvKfuaux7an+vFSFgo5HdHhJ9rYi+uH+VovPXO+64PBo56j/R6y/8n6prvXfM/s5ErQKNbebeX7lq1n2/lsZVNowgdFsiqlnrTtQXfFzccX0ae/4Syk+qE8jkG1FzVS5R5D2YqCbiZKJRSoTt84Fz0pjWWUTV+A26eA4gndPDB1i7FJapPdEItnSevzOIfrk3pi9LIuWaD5yXpjKdnYbz156oQ/cZaXytGnVdEDXLtiWujx7RF8BvgP8K9/OTUKZE02r0L/YDTauh1/oBpHNY+NvgoKqZnr/Dw98G9WU+gHQmpbFM4zM0rYlJSKcNMX1R2RsUPgQMTdZrdSA39YFrIIsmqXyBqA/TX4n6nT0KXA18ysz6unspUbVug/srNCAt3L3C3Q+4P9ABppPl0eiv0hSWaTdRn4Z0nj/zSHmKy/XPsC8dZTrgPk0NSGs30fx2jZKu91VI53mipt97zWx0zT53f5Jo/saeZjYTeMkaMZ3MAab1ak3f2TSk9e9wPg/4Wj/AdN4I/TEb9Ho14Pyl67V6La6faarSeTmNZfp3hqb1SkPfVzHp/J5oQuCadGr6K3cDOpjZZcA0S+50RvVLZ7SYCTei2pmBRO3bJxA1q1xPNEKyP1ENzgPhdgewDhjZnNPKxDJlalqZWKZmcP6+TzQ339i4Yx8i6ofW0P4yGZdWJpZJ509pHUA618anQ/RD8wmieeYa1W3pgPOYzsQy5UbU9n030RdLTTXqNUSjPfsS9fc4h2hUV6NWREhXWplYpkxNKxPL1AzO39VEoyVHhvt9gXdoZDNPJqaViWXS+VNaDUhnVLj/K2AVDRwE1KjzkO4EW/KNaKjzEUQdFx8jmgAxdv+NRNF4o+Y4SmdamVimTE0rE8vUzM7f9URrT9b0tWvwckKZmFYmlknnT2k1Ip37iQK804BhDS1TY25pT7Cl3ojWeFxC1K7+W6JJS1cDN8YcM4QoWm/wBIHpTCsTy5SpaWVimZrp+ftTGsvUYtLKxDLp/CmtRqbz58aUJxm3Jk28pdyAo4mW66kZtXU30WoA/YjmdfsRUcR+JdEcMY2ZTT8taWVimTI1rUwsk85fy0krE8uk86e0kpROo1aOaOytyRJuSbfwol4Zc78ne5fwGUbUvPN7GjlJZTrTysQyZWpamVgmnb+Wk1YmlknnT2k1dZmScWvSxFvKjaidu3PM/wOI5rzqG7YNJhpS3KWlpJWJZcrUtDKxTDp/LSetTCyTzp/SauoyJeOmeeAS4O5V7r493DWgGNjm7uvD3C8/JJrgr8FrMaY7rUwsU6amlYllSmdamVimdKaViWVKZ1qZWKZMTSudZUqGmmGxcoDM7D6i+WBOIapyXdrS08rEMmVqWplYpnSmlYllSmdamVimdKaViWXK1LTSWaYDpQDuAJmZES2psSL8/Yy7/6clp5WJZcrUtDKxTOlMKxPLlM60MrFM6UwrE8uUqWmls0wNpQCugczsSmCuuy/PlLQysUyZmlYmlimdaWVimdKZViaWKZ1pZWKZMjWtdJbpQCmAayAzM0/TyUtXWplYpkxNKxPLlM60MrFM6UwrE8uUzrQysUyZmlY6y3SgFMCJiIiItDAahSoiIiLSwiiAExEREWlhFMCJiIiItDAK4ERERERaGAVwIiIiIi2MAjgRkTQys+ymzoOItHwK4ERE6mBmt5jZd2Pu325mV5vZdWY218yWmNlPY/bPNLP5ZrbczK6K2b7TzH5lZouBo9JbChHJRArgRETq9lfgCgAzywIuBjYAI4AjgfHAJDM7Lhz/JXefBOQD3zGz7mF7B+Atdz/c3V9LY/5FJEPlNHUGRESaK3dfbWZbzWwC0BtYCBxBtLD1wnBYR6KA7lWioO28sH1g2L4VqAKeTGfeRSSzKYATEanfX4ArgT5ENXKfAe5w9z/FHmRmJwAnAUe5+24zexnIC7vL3L0qTfkVkVZATagiIvWbAUwlqnmbHW5fMrOOAGbW38x6AV2Aj0LwNhqY0lQZFpHMpxo4EZF6uPseM3sJKA61aH8zs0OAN80MYCdwGfAi8HUzWwG8C8xpqjyLSObTYvYiIvUIgxcWAJ919/80dX5EREBNqCIidTKzMcAq4J8K3kSkOVENnIiIiEgLoxo4ERERkRZGAZyIiIhIC6MATkRERKSFUQAnIiIi0sIogBMRERFpYf4/4cNrcnrSDJsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 可视化\n",
    "fig, ax = plt.subplots(figsize=(10, 4))\n",
    "ax.plot(data.index, data['value'], 'o-')\n",
    "ax.set_title(u\"Data on urban employees in the past 20 years\")\n",
    "ax.set_ylabel(\"Ten thousand people\")\n",
    "ax.set_xlabel(\"year\")\n",
    "ax.set_xticks(data.index)\n",
    "ax.set_xticklabels(data.index, rotation=45)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# lag（滞后）\n",
    "经常在时间序列里面可以看到lag这个东西，大白话就是叫滞后；\n",
    "举个例子来说：\n",
    "\n",
    "1. lag1代表滞后数据1项，也就是2002年的数据滞后1年得到的就是2001年的数据；2003年的数据滞后1年就是2002年的数据；\n",
    "2. lag2代表滞后数据2项，也就是2003年的数据滞后2年得到2001年的数据；2004年的数据滞后2年就是2002年的数据；\n",
    "3. 依次类推，滞后数据就是这么计算的。\n",
    "\n",
    "下面使用python对上面的【近20年城镇就业人员数据(万人)】数据做个滞后处理。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>value_lag1</th>\n",
       "      <th>value_lag2</th>\n",
       "      <th>value_lag3</th>\n",
       "      <th>value_lag4</th>\n",
       "      <th>value_lag5</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>year</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2001</th>\n",
       "      <td>24123</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2002</th>\n",
       "      <td>25159</td>\n",
       "      <td>24123.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2003</th>\n",
       "      <td>26230</td>\n",
       "      <td>25159.0</td>\n",
       "      <td>24123.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004</th>\n",
       "      <td>27293</td>\n",
       "      <td>26230.0</td>\n",
       "      <td>25159.0</td>\n",
       "      <td>24123.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2005</th>\n",
       "      <td>28389</td>\n",
       "      <td>27293.0</td>\n",
       "      <td>26230.0</td>\n",
       "      <td>25159.0</td>\n",
       "      <td>24123.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2006</th>\n",
       "      <td>29630</td>\n",
       "      <td>28389.0</td>\n",
       "      <td>27293.0</td>\n",
       "      <td>26230.0</td>\n",
       "      <td>25159.0</td>\n",
       "      <td>24123.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2007</th>\n",
       "      <td>30953</td>\n",
       "      <td>29630.0</td>\n",
       "      <td>28389.0</td>\n",
       "      <td>27293.0</td>\n",
       "      <td>26230.0</td>\n",
       "      <td>25159.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2008</th>\n",
       "      <td>32103</td>\n",
       "      <td>30953.0</td>\n",
       "      <td>29630.0</td>\n",
       "      <td>28389.0</td>\n",
       "      <td>27293.0</td>\n",
       "      <td>26230.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2009</th>\n",
       "      <td>33322</td>\n",
       "      <td>32103.0</td>\n",
       "      <td>30953.0</td>\n",
       "      <td>29630.0</td>\n",
       "      <td>28389.0</td>\n",
       "      <td>27293.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2010</th>\n",
       "      <td>34687</td>\n",
       "      <td>33322.0</td>\n",
       "      <td>32103.0</td>\n",
       "      <td>30953.0</td>\n",
       "      <td>29630.0</td>\n",
       "      <td>28389.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      value  value_lag1  value_lag2  value_lag3  value_lag4  value_lag5\n",
       "year                                                                   \n",
       "2001  24123         NaN         NaN         NaN         NaN         NaN\n",
       "2002  25159     24123.0         NaN         NaN         NaN         NaN\n",
       "2003  26230     25159.0     24123.0         NaN         NaN         NaN\n",
       "2004  27293     26230.0     25159.0     24123.0         NaN         NaN\n",
       "2005  28389     27293.0     26230.0     25159.0     24123.0         NaN\n",
       "2006  29630     28389.0     27293.0     26230.0     25159.0     24123.0\n",
       "2007  30953     29630.0     28389.0     27293.0     26230.0     25159.0\n",
       "2008  32103     30953.0     29630.0     28389.0     27293.0     26230.0\n",
       "2009  33322     32103.0     30953.0     29630.0     28389.0     27293.0\n",
       "2010  34687     33322.0     32103.0     30953.0     29630.0     28389.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 使用pandas的shift可以做lag\n",
    "data1 = data.copy()\n",
    "# data1['value_lag1'] = data1['value'].shift(1)\n",
    "# data1['value_lag2'] =  data1['value'].shift(2)\n",
    "for i in [1, 2, 3, 4, 5]:\n",
    "    data1[f\"value_lag{i}\"] = data1['value'].shift(i)\n",
    "data1.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "1. 因为2001年之前是没有数据了，所以2001年滞后一项得到的数据是空数据；\n",
    "2. 因为2002年之前是没有数据，所以2002年滞后两项也是没有数据的。\n",
    "3. 依次类推。\n",
    "\n",
    "# 时间序列的公式\n",
    "\n",
    "经常看到时间序列的公式写的乱七八糟的，这是一堆公式，那是一堆公式，活活把人绕晕。现在为了我们说明需要，我们定义几个简单的公式：\n",
    "\n",
    "假设我们的时间序列是`value = [4,3,2,10,3]`。我们把这个时间序列定义为 $ X_{t} $；\n",
    "那么这个时间序列里面每一个元素对应的符号依次是：$x_1$、$x_2$、$x_3$等。 \n",
    "\n",
    "如果$X_t$做了lag1，得到的时间序列为`value lag1=[nan, 4,3,2,10]`，记这个新的时间序列为$X_{t-1}$。\n",
    "\n",
    "如果$X_t$做了lag2，得到的时间序列为`value lag2=[nan,nan, 4,3,2]`，记这个新的时间序列为$X_{t-2}$。\n",
    "\n",
    "以此类推，就是一个简单的时间序列公式了。\n",
    "\n",
    "# 统计基础知识\n",
    "\n",
    "1. 均值:$\\mu = \\frac{\\sum X}{n}$\n",
    "\n",
    "\n",
    "2. 方差：$S^2 = \\frac{\\sum(X-\\mu)^2}{n-1}$\n",
    "\n",
    "\n",
    "3. 协方差：$Cov(X,Y) = E\\{|EX - X| |EY-Y|\\}$\n",
    "\n",
    "\n",
    "4. 相关系数: $\\gamma(X, Y) = \\frac{Cov(X, Y)}{\\sqrt{Var(X)Var(Y)}}$\n",
    "\n",
    "# 自相关系数\n",
    "\n",
    "## 自相关基础计算\n",
    "在了解上面的【统计基础知识】后，再理解自相关性基本上就不难了，自相关是变量与自身在不同时间滞后下的相关性。对于延迟$k$阶的时间序列。自协方差为$c_k$，自相关性系数为$\\gamma_k$\n",
    "\n",
    "1. $c_k = \\frac{1}{n}\\sum_{i=1}^{n-k}(x_t - \\overline{x})(x_{t+k}-\\overline{x})$\n",
    "\n",
    "\n",
    "2. $\\gamma_k = \\frac{c_k}{c_0}$\n",
    "\n",
    "### 调用包计算acf数值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 使用statsmodels包的acf函数\n",
    "from statsmodels.tsa.stattools import acf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 自己根据公式写函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cal_acf(x, nlags):\n",
    "    \"\"\"\n",
    "    按照公式自己写个acf函数\n",
    "    只是用来帮助理解，不建议用在别的环境中\n",
    "    \"\"\"\n",
    "    x = np.array(x)\n",
    "    mean_x = np.mean(x)\n",
    "    length_x = x.shape[0]\n",
    "    c_0 = np.mean((x-mean_x) **2)\n",
    "    c_k = np.sum((x[:(length_x-nlags)] - mean_x) * (x[nlags:] - mean_x)) / length_x\n",
    "    r_k = c_k / c_0\n",
    "    return r_k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 两种结果对比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/devc/lib/python3.8/site-packages/statsmodels/tsa/stattools.py:667: FutureWarning: fft=True will become the default after the release of the 0.12 release of statsmodels. To suppress this warning, explicitly set fft=False.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>value_by_myself</th>\n",
       "      <th>value_by_statsmodels</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.858910</td>\n",
       "      <td>0.858910</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>0.715601</td>\n",
       "      <td>0.715601</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>0.570864</td>\n",
       "      <td>0.570864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>0.427266</td>\n",
       "      <td>0.427266</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>0.287384</td>\n",
       "      <td>0.287384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "      <td>0.154147</td>\n",
       "      <td>0.154147</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>7</td>\n",
       "      <td>0.031085</td>\n",
       "      <td>0.031085</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>8</td>\n",
       "      <td>-0.082773</td>\n",
       "      <td>-0.082773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>9</td>\n",
       "      <td>-0.184306</td>\n",
       "      <td>-0.184306</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>10</td>\n",
       "      <td>-0.269526</td>\n",
       "      <td>-0.269526</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    index  value_by_myself  value_by_statsmodels\n",
       "0       0         1.000000              1.000000\n",
       "1       1         0.858910              0.858910\n",
       "2       2         0.715601              0.715601\n",
       "3       3         0.570864              0.570864\n",
       "4       4         0.427266              0.427266\n",
       "5       5         0.287384              0.287384\n",
       "6       6         0.154147              0.154147\n",
       "7       7         0.031085              0.031085\n",
       "8       8        -0.082773             -0.082773\n",
       "9       9        -0.184306             -0.184306\n",
       "10     10        -0.269526             -0.269526"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 结果汇总\n",
    "pd.DataFrame({'index':np.arange(11),\n",
    "             'value_by_myself':[cal_acf(x=data['value'], nlags=i) for i in range(11)],\n",
    "             'value_by_statsmodels':acf(data.value,nlags=10)})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 自相关系数画图\n",
    "\n",
    "画图有两种方法：\n",
    "1. 直接使用`statsmodels`的`plot_acf`;\n",
    "2. 使用`acf`函数来计算数据，得到数据后再进行可视化，这种一般便于自己定制，也方便自己做更多的细节调整之类的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0, 15.5)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0YAAAIHCAYAAAC7V0xlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVpklEQVR4nO3deZhcV33n//e3W6utxYvakvBuwMIBE7NvQbYBzzhhJgk7CUkgmUzMYoKBEDAkMYZMDCEYMzgEDGHLkLAa50cS8ASCkAGPiQFjsQnjBWRbsiXZ6tbSa/X5/XFuqUul3tVVt5b363nuU1237q36drWW+vQ553sjpYQkSZIkdbOesguQJEmSpLIZjCRJkiR1PYORJEmSpK5nMJIkSZLU9QxGkiRJkrqewUiSJElS1zMYSZIkSep6BiNJkiRJXc9gJEmSJKnrGYwkSQBExNqI+FxE7I6IFBGXLMBzvjUi0gKUN9PrnFbU/CeNfq12ExHnFe/NeXX7fzcifhIRoxGxp5TiJKmFLCq7AElSy3gP8F+By4EdwM1lFhMRrwQOpJQ+VmYdnSgiHgF8DPgy8A7gQKkFSVILMBhJkqqeAfxzSulvyi6k8EpgF/kDvBbWeeRZI69JKf2s5FokqSU4lU6SVHUCsKfsItQUJxS3e8osQpJaicFIktpURJwaEe+PiK0RMVisDfpsRJw2ybHHRMR7IuKuiBiOiLsj4hMRsSYiXlasAwrgVcV6lCnXBdWu54mI10bEz4vX/3pEPGoWdS+KiD+PiNuLWu6KiL+KiKU1x9wFPBI4t1pPRGya5fsyZU0R8fvFcz1mkvPeHBGViDhxmud+a3H+mRHxfyKiPyJ2RsTbIzs5Iv45IgYiYkdEvL7m3BURsT8i3jvJ855UvPalxf3FEXFZRNwWEUPFz/YbEXFB3XmPKNaFPVAcd3NE/PoM789d5OmSADuL7+et050jSd3AqXSS1L6eADwV+BRwN3Aa8ApgU0T8UkrpAOQP5MANwFnAR4DvAmuAXwdOAjYDvwv8A/DvwCdm+fq/B6wE/hZYBrwG+I+IODuldN80530YeCnwOeDdwJOAS4v6nlMccwnwPmAf8L+KfdM952xr+lzx2EuA79Wd+xJgU0rpnlm8zqeBHwNvAp4N/BnwAHAR8B/AG4vn+5uI+M+U0uaU0r6I+ALwooh4XUqpUvN8v0UOpp8s7r+V/J58GPg2sAp4PPBY8s+IiHgk8E3gHvI6of3AC4HrIuJ5KaUvTFH7JcX79Bzyn5d9wK2z+J4lqbOllNzc3Nzc2nADlk+y78lAAn63Zt/lxb7nTHJ81HydgKtn8bqnFcceAE6s2f/EYv+VNfvemv+rOXj/l4tjPlT3nO8q9p9fs+8H5KAym/diLjX9IzlM9NTse0xx3MtmeJ23Fsd9sGZfL7ANGAfeWLP/mKKej9Xs+y/F+RfWPe/3a79X4BbgX2ao5SvkQLO09udJDks/rdl3XvGa503yfawp+8+xm5ubW6tsTqWTpDaVUhqsfl1MvToe+Bl53chjaw59HvD9NMkIQkrpSFppX5dqRldSSt8GbgJ+bZpzqo9dWbf/3cXts4+gntnW9AngIcD5NfteAgwCn5/l63y45jUq5A5+Afx9zf49wFbgjJrzvgLcW7weAMVUv0cD/6fmuD3AIyPi4ZO9eEQcR26W8RlgZTElcg1wPHA98PDppgRKkg5nMJKkNhURyyPibRGxDRgmd3DbSR6pWF1z6EPJoy8L7bZJ9v2UPHozlVPJIyuHdEJLKe0gh4FTm1DTvwPbKcJJRPSQp7L9c0pp7yxf5xd19/uBoZTSrkn2H1u9k1IaJ0+X+82IOKrY/RJgCPhszXl/Qf45/jQitkTEuyLi0TWPP4wcxN5O/pnXbtX1QycgSZo1g5Ekta/3AW8hjxq8kDxN6wJgN63/73vDL/o65QvnEZ5/BJ4XEcvII0cP4dARm5lUZrkPcoCp9QlgBTkcBfDb5Glz/TU1biYH2j8gh9o/BL4bEX9YHFL9+f4N+Wc+2WYbbkmaA5svSFL7ej7w8ZRSbeezZeSRhlq3AzN2i5uHyaZ5nQncNc05Pyd/qH84uXkBABGxllz3z2uOnU94mm1NnwBeD/x34FfJIy3Xz+P15iyl9IOI+B55pOhu4BTg1ZMc9wDwUeCjRQONzeS1QR8G7igOG00pfaUZdUtSp2v13yhKkqZW4fDRiFeTmwHU+jzwyxHxnLr9FCMW8/WbtetYIuKJ5A5zX5rmnH8rbi+p2/+64vZfa/bt5/CQtyA1pZRuJTcu+EPyGqxPpZTG5vhaR+IfyCN8l5BH+A6pr1gvdlBKaR95BGhpcf9+YBNwUUSsr3/yiOhrRNGS1MkcMZKk9vUvwO9GRD/wI+ApwLPIH7RrvYs8uvTZiPgI8B3gOHK77peTO6LNx8+Ab0TE35E/sF9SvPZfT3VCSun7EfFx4I8i4hjg6+TOcS8lN074Ws3h3wFeERF/VrzW/Sml/1jAmj5BnooGc5tGtxD+kVzTc4C/SymN1j3+o+K6Td8htwF/PPlneHXNMa8CvgFsiYgPkUeR1pL/HJxE7gAoSZolg5Ekta/XkEeNXkK+Zs83ycHokClhKV8/5+nkRfnPIYeQ+4GvkqdyzdcnyI0ULiEv9P82cHFKafsM5/0h+UP8y4p6dgBXMNE0oOpt5GYMf0q+NtHXydcIWqiaPgm8E7i96F7XNCml+yLi/5K75f3DJIf8b3Jw/S/kgPdz8rWS3lXzHD+KiMcDl5Hfy+PJP9fvkd87SdIcxJF1apUkdZuIOA24E3hDSulvZji8ZRXtrbcDb0spvb2E1/8CcHZK6WHNfm1J0uFcYyRJ6lYvI6/HmmzEpqGKdUHPLuO1JUmTcyqdJKmrRMQzgF8itzq/LqV0VxNf+3TgaeTphKPAB5v12pKk6RmMJEnd5i+Ap5LXZB3WJrvBziW34P4F8NLiwraSpBbgGiNJkiRJXc81RpIkSZK6nsFIkiRJUtfruDVGxVXcHwLsLbsWSZIkSaVbCdybZlhD1HHBiByKjuSChZIkSZI6y0nAPdMd0InBaC/Atm3bWLVqVdm1SJIkSSrJwMAAJ598MsxiNlknBiMAVq1aZTCSJEmSNCs2X5AkSZLU9QxGkiRJkrqewUiSJElS1zMYSZIkSep6BiNJkiRJXc9gJEmSJKnrGYwkSZIkdb2GXscoIjYCbwAeB6wHnpNSum6Gc84DrgQeCWwD/jKl9LFG1tkSKhW44QbYvh3Wr4enPx16e8uuSpIkSe0mpfzZslKB8fGFuV20CM45p+zvrKEafYHXo4HvAx8Brp3p4Ig4HfhX4APAS4BnAh+OiO0ppesbWWiprr0WXvMauPvuiX0nnQTvfS8897nl1SVJkqSFNVMAWYgQk9LC171s2cI/Z4tpaDBKKX0J+BJARMzmlJcDd6aUXl/c/3FE/ArwWqAzg9G118Lzn3/4H+B77sn7P/c5w5EkSVIz1IaTuQSR2R47Pt6Y0KIF0egRo7l6CvCVun3XA1c1v5QmqFTySNFkf0FSggi45BL4jd9wWp0kSdJkqsFjbOzQ26m+nm6foaWrtVowWgfcV7fvPmBVRCxPKQ3WnxARS4GlNbtWNrC+hXXDDYdOn6uXEmzblo8777ymlSVJktRQtSMtswks0309Pl72d6MO0WrBaD4uBS4ru4h52b59YY+TJElaSJMt4p/t19OFHMOMWlCrBaMdwNq6fWuBgclGiwpXkLvYVa0EphmGaSHr1y/scZIkqbscSXCZzddOLVMXabVgdCPwa3X7Lij2TyqlNAwMV+/PsslDa3j603P3uXvumfwfnoj8+NOf3vzaJEnSzFKaWFRfXVg/2/uzOdbgIjVNo69jtAJ4WM2u0yPiHOCBlNIvIuIK4MSU0u8Vj38AuDgi/prc4vsZwAuBZzeyztL09uaW3M9/fg5Btf+4VQPeVVfZeEGS1DlSmvj/rvr1dPsW4tiFDiu19yV1jEaPGD0e+FrN/eqUt48DLyNf9PWU6oMppTsj4tnAe4DXkKfE/WFHX8Pouc/NLbknu47RVVfZqlvSkaltEVv7wW6yD5LV28n2zfexRhw/G636W/Sp6mqH/QsVYCSpRUXqsH+kImIV0N/f38+qVavKLmf2KpXcfe71r4d3vztPn3OkSGp/1akw9eGkWZskSQth2TJ48pPLrmLOBgYGWL16NcDqlNLAdMe22hqj7tXbm1tyn3iirbmlslUqMDqaOyjV3tbvm03g6bBfPkmS1KkMRjpy1dGu7dtzBz1Hu9QqxsenDzb1+6pfO9IiSVLXMRjpyFx77eTro977XtdHaeGkNHOYmWpER5IkaRYMRpq/a6/NHfXqpwrdc0/e/7nPGY40uZRgaAgGBycCzXRhx4AjSZIazGCk+alU8kjRVJ2LIuCSS+A3fsNpdd2qUpkIP/Xb8LBrbyRJUksxGGl+brjh0Olz9VKCbdvycTaT6FyjoxNhpz4EjYyUXZ0kSdKsGYw0P9u3L+xxal3Dw4cGntoANDZWdnWSJEkLwmCk+Vm/fmGPU3lSmnzEp7rPDm2SJKkLGIw0P09/eu4+d889k68ViciPP/3pza9Nh6tUpp7y5nofSZIkg5Hmqbc3t+R+/vNzCKr9YB2Rb6+6ysYLzVS73qc+BLneR5IkaVoGI83fc5+bW3JPdh2jq66yVXcjjY1Bf3/e9uyBAwdc7yNJknQEDEY6Ms99bm7JfcMN8PrXw7vfnafPOVK0sEZHJ0JQfz/s2+f0N0mSpAVkMNKR6+3NLblPPNHW3AtldDSHoNogJEmSpIYxGEmtYGRkIghVp8ZJkiSpaQxG6m6VSp4GuH17bi3erGmAw8OHBqHBwca/piRJkqZkMFL3uvbayRtHvPe9C984YnBwYo3Qnj25W5wkSZJahsFI3enaa3Or8foGBvfck/d/7nNHFo4OHDg0CA0PH0m1kiRJajCDkbpPpZJHiibr6pZSvg7TJZfkbnuznVa3f/9Eo4Q9e7xukCRJUpsxGKn73HDDodPn6qUE27bl4ybrspfS4UFodLRBxUqSJKkZDEbqPtu3z+24lHK77Nr22V5MVZIkqaMYjNR91q+f3XE9PXDrrTAwYBCSJEnqcD1lFyA13dOfnrvPRUx9TF8frFkDDzxgKJIkSeoCBiN1n95euOqqyZsvVF18cXOuZyRJkqSWYDBS96hUYOdO+PGPYe1auPzyPCpUq68v79+4sZwaJUmSVArXGKmzjY7C7t05ED34IIyPTzy2cSM87WmwZQu8//3wylfC2Wc7UiRJktSFDEbqPENDsGtX3vr7p58y19sL55yTR4rOOadZFUqSJKnFGIzUGfbvnwhDe/eWXY0kSZLajMFI7WtgIE+R27ULBgfLrqb5KpU8DXD3bjj+eKcBSpIkHQGDkdpHSnmdUHVkaGSk7IrKs3kzvO99+X2oWrMGXv1qG0dIkiTNg8FIra1SydcS2rUrj4x4TaEcii677PD9u3bl/XbVkyRJmjODkVrP6OjEqFB9J7luV6nkkaLpXH117rbntDpJkqRZMxipNcylk1w327Ll0Olzk9m5Mx9nlz1JkqRZMxipPNVOcjt3wr59ZVfTHnbvXtjjJEmSBBiM1Gz9/RMjQ93YSe5IHX/8wh4nSZIkAHqa8SIR8aqIuCsihiLipoh44jTHviwiUt021Iw61QAp5eYJP/0pfOtb8L3vwbZthqL5Ovvs3H1uOn19+ThJkiTNWsNHjCLiRcCVwMuBm4BLgOsjYkNK6f4pThsANtTcd8FJO6l2ktu5M9/aSW7h9PbmltyTdaWruvhiGy9IkiTNUTNGjF4HfCil9NGU0o/IAekA8AfTnJNSSjtqtvuaUKeOxOhobqCwZQt885vwwx/C/fcbihph48bckrt+5Kivz1bdkiRJ89TQEaOIWAI8Driiui+lNB4RXwGeMs2pKyLi5+Tg9l3gzSmlH07xGkuBpTW7Vh5x4Zqd+k5y+/e76L9ZNm7MLbm3bIH3vx9e+co8fc6RIkmSpHlp9FS6NUAvUD/icx/wiCnO2UoeTboVWA38CfCtiHhkSunuSY6/FJhmXpEWlJ3kWkdvb27J3ddna25JkqQj1HJd6VJKNwI3Vu9HxLeAHwMXAX8+ySlXkNcwVa0EJgtQmi87yWmhVSp5tGv37txBz9EuSZJUskYHo11ABVhbt38tsGM2T5BSGo2I7wEPm+LxYWC4ej8i5lepJqQEDz44EYZGRsquSJ1k82Z43/sOvVDtmjW5qYTroyRJUkka2nwhpTQCfAd4ZnVfRPQU92+c6rxaEdELnA1sb0SNKlQqeXrcj36Umyfceivce6+hSAtr8+bcUa82FEG+f9ll+XFJkqQSNGMq3ZXAxyPiZuDb5HbdRwMfBYiITwD3pJQuLe7/BfD/gJ8BxwBvAE4FPtyEWrvL6OjEqNCDD8L4eNkVqZNVKnmkaDpXX52bSjitTpIkNVnDg1FK6dMR0Qe8DVgH3AJcWNOC+xSg9hP5scCHimMfJI84PbVo9a0jVd9JLnmJKDXJli2HjxTV27kzH2czCUmS1GRNab6QUroauHqKx86ru/9a4LVNKKs73XabLbVVjtn+ufPPpyRJKkEzLvAqSbn73EIeJ0mStIAMRpKa4+yzc/e56fT15eMkSZKazGAkqTl6e3NL7ulcfLGNFyRJUikMRpKaZ+NGuPzyw0eO+vryfq9jJEmSStKU5guSdNDGjbkl95Yt8P73wytfmafPOVIkSZJKZDCS1Hy9vbkld1+frbklSVJLcCqdJEmSpK7niJEkzUWlkqcB7t6dW4s7DVCSpI5gMJKk2dq8Gd73Pti1a2LfmjW5256NIyRJamtOpZOk2di8GS677NBQBPn+ZZflxyVJUtsyGEnSTCqVPFI0nauvzsdJkqS2ZDCSpJls2XL4SFG9nTvzcZIkqS0ZjCRpJrt3L+xxkiSp5RiMJGkmxx+/sMdJkqSWYzCSpJmcfXbuPjedvr58nCRJaksGI0maSW9vbsk9nYsv9npGkiS1MYORJM3Gxo1w+eWHjxz19eX9XsdIkqS25gVeJWm2Nm6Epz0td597//vhla/M0+ccKZIkqe0ZjCRpLnp74Zxz8kjROeeUXY0kSVogBiNJ6haVSh7t2r07d9BztEuSpIMMRpLUDTZvhve979AL1a5Zk5tKuD5KkiSbL0hSx9u8GS677NBQBPn+ZZflxyVJ6nIGI0nqZJVKHimaztVX5+MkSepiBiNJ6mRbthw+UlRv5858nCRJXcxgJEmdbPfuhT1OkqQOZTCSpE52/PELe5wkSR3KYCRJnezss3P3uen09eXjJEnqYgYjSepkvb25Jfd0Lr7Y6xlJkrqewUiSOt3GjXD55YePHPX15f1ex0iSJC/wKkldYeNGeNrTcve5978fXvnKPH2uXUaKKpVc++7deT1UO9UuSWoLBiNJ6ha9vXDOOXmk6Jxzyq5m9jZvztdiqm07vmZNniLoaJckaYE4lU6S1Lo2b4bLLjv8Wky7duX9mzeXU5ckqeMYjCRJralSySNF07n66nycJElHqCnBKCJeFRF3RcRQRNwUEU+c4fgXRMRPiuO3RMSvNaNOSVIL2bLl8JGiejt35uMkSTpCDQ9GEfEi4ErgcuCxwPeB6yPihCmOfyrwT8DfA48BrgOui4hHNbpWSVIL2b17YY+TJGkazRgxeh3woZTSR1NKPwJeDhwA/mCK418DfDml9K6U0o9TSn8OfBe4uAm1SpJaxfHHL+xxkiRNo6HBKCKWAI8DvlLdl1IaL+4/ZYrTnlJ7fOH6qY6PiKURsaq6ASuPuHBJUvnOPvvway/V6+vLx0mSdIQipdS4J494CHAP8NSU0o01+/8aODel9KRJzhkBXppS+qeafa8ELksprZ3k+LcCl9Xv77/wQlYtXrwg30dT3XMPnHhi455/YABGRxvz3Dt35g8p7cjay2Ht5Win2nfvhq1bp358w4b2GTFqp/e9nrWXw9rLYe2T6+mBY49tzHM30MDoKKu//GWA1SmlgemO7YTrGF1BXsNUtRK4m09/GlatKqmkFla9QKIktYvJrmPU1wcXX9xe1zF6y1vgf/2vsquYH2svh7WXw9ont2wZPPnJjXnuRhoYgNWrZ3Voo4PRLqAC1I/0rAV2THHOjrkcn1IaBoar9yNiXoVKklrUxo3wtKdN/GLn+OPz9Lne3rIrkyR1kIauMUopjQDfAZ5Z3RcRPcX9G6c47cba4wsXTHO8JKnT9fbCOefAM5+Zb9spFFUqcMsteYrLLbd43SVJalHNmEp3JfDxiLgZ+DZwCXA08FGAiPgEcE9K6dLi+PcCX4+I1wP/CrwYeDzwR02oVZKkhVM/DfC1r80NJV796vaaBihJXaDh7bpTSp8G/gR4G3ALcA5wYUrpvuKQU4D1Ncd/C/htchD6PvB84DdTSj9odK2SJC2YzZvhsssOv0jtrl15/+bN5dQlSZpUU5ovpJSuBq6e4rHzJtn3WeCzDS5LkqTGqFTySNF0rr46r51qp2mBktTBmnGBV0mSusuWLYePFNXbuTMfJ0lqCQYjSZIW2mwvi+DlEySpZRiMJElaaLO96Gy7XJxWkrqAwUiSpIV29tm5+9x0+vrycZKklmAwkiRpofX25pbc07n4YhsvSFILMRhJktQIGzfC5ZcfPnLU15f3t8N1jLw4raQu0pR23ZIkdaWNG3NL7i1bcqOF44/P0+faYaTIi9NK6jIGI0mSGqm3F845p+wq5qZ6cdp61YvTtsuIlyTNgVPpJEnShNlenNZpdZI6jMFIkiRN8OK0krqUwUiSJE3w4rSSupTBSJIkTfDitJK6lMFIkiRN6JSL09pqXNIcGYwkSdKETrg47ebN8OIX5xbjt92Wb1/84rxfkqZgMJIkSYdq54vTVluN1zeQqLYaNxxJmoLXMZIkSYdrx4vTzrbV+NOe1trfh6RSGIwkSdLk2u3itHNpNd5O35ekpnAqnSRJ6gy2Gpd0BAxGkiSpM3RKq3E76kmlcCqdJEnqDNVW49NNp2v1VuObN+d1UtXv4bWvzd/Tq1/d2k0vpA7giJEkSeoM7d5q3I56UqkMRpIkqXO0a6vx2XbUc1qd1DBOpZMkSZ2lHVuNd0pHvUol11hdH9Xq77tUw2AkSZI6T7u1Gu+Ejnquj1KbcyqdJElS2dq9o57ro9QBDEaSJEllq3bUm06rdtTrlPVRtknvegYjSZKksrVzR725rI9qVZs3w4tfnKf/3XZbvn3xi9tnpMtQtyAMRpIkSa2gXTvqtfv6qHafBtjuoa6F2HxBkiSpVbRjR712Xh8122mAT3taa/4MqqGuXjXUtXKgbkGOGEmSJLWSake9Zz4z37biB/Ja7bw+qp2nAXbK2q4WYjCSJEnS/LXz+qh2ngbYzqGuRRmMJEmSdGTadX1UO08DbOdQ16JcY9Rtjj0W9u2D4eGyK5EkSZ2kHddHVacBTjfy0qrTANs51LWoho4YRcRxEfHJiBiIiD0R8fcRsWKGczZFRKrbPtDIOrvKSSfBU54Cj30snHIKHHVU2RVJkqRO0W7ro9p5GmA7r+1qUY2eSvdJ4JHABcB/AzYC18zivA8B62u2P21UgV1r1So44wx44hPzdvrpsHJl2VVJkiQ1V7tOA2znUNeiIqXUmCeOOAv4EfCElNLNxb4LgX8DTkop3TvFeZuAW1JKl8zzdVcB/f39/axatWo+T9HdhofzcPKuXbBnDzToz4ckSVJLqVTaaxpg1ebNuTtd7XTAvr4cihYy1C1bBk9+8sI9X5MMDAywevVqgNUppYHpjm1kMPoD4N0ppWNr9i0ChoAXpJS+MMV5m8ijTAHsAL4IvD2ldGCK45cCS2t2rQTuNhgtgNHR/I/Drl3wwAMwPl52RZIkSarXjFDXBcGokc0X1gH31+5IKY1FxAPFY1P5R+DnwL3Ao4F3AhuA505x/KXAJFe20hFbvBjWrctbpQIPPpjbPu7eDWNjZVcnSZIkmFjbpSMy52AUEe8A3jjDYWfNrxxIKdWuQdoSEduBr0bEQ1NKt09yyhXAlTX3VwJ3z/f1NYXe3jz3ds2aPL1uz56JKXd2uJMkSVKbm8+I0buBj81wzB3kaXAn1O4sptIdVzw2WzcVtw8DDgtGKaVh4OAn84iYw1NrXiJy2+9jj4WHPxwGBiZC0oFJZzxKkiRJLW3OwSiltBPYOdNxEXEjcExEPC6l9J1i9zPInfBumvrMw5xT3G6fS51qolWrJrrcHTiQp9vt2gV795ZdmSRJkjQrDVtjlFL6cUR8GfhQRLwcWAxcDXyq2pEuIk4Evgr8Xkrp2xHxUOC3yZ3rdpPXGL0H2JxSurVRtWoBHXUUnHpq3uxwJ0mSpDbRyOYLAC8hh6GvAuPA54E/rnl8MbmxQvUqoyPAs4BLgKOBbcU5f9ngOtUIS5fCiSfmzQ53kiRJamENDUYppQfII0BTPX4XuS139f424NxG1qSS2OFOkiRJLazRI0bS4exwJ0mSpBZjMFK57HAnSZKkFmAwUmup73C3a1eecmeHO0mSJDWQwUit66ij4JRT8maHO0mSJDWQwUjtwQ53kiRJaiCDkdrPZB3uqqNJdriTJEnSPBiM1N7scCdJkqQFYDBS56jvcLd3b27cYIc7SZIkzcBgpM61cmXeajvc7dqVW4JLkiRJNQxG6g52uJMkSdI0DEbqPrUd7sbGcoe7nTtzSLJ5gyRJUlcyGKm7LVoEa9fmLSXYvz8HpD17oL8/twaXJElSxzMYSVURsGJF3k46Ke+rD0ojI2VWKEmSpAYxGEnTOfrovJ14Yr5/4MBESNqzx5bgkiRJHcJgJM3FUUfl7SEPyfcHBw8NSkNDZVYnSZKkeTIYSUdi+fK8rV+f7w8NHRqUBgfLrE6SJEmzZDCSFtKyZbBuXd4gT7WrhqQ9e7zQrCRJUosyGEmNtHQpnHBC3iA3b6gNSvv3l1mdJEmSCgYjqZmWLIG+vrxBbgdeH5S84KwkSVLTGYykMi1eDGvW5A3yBWarQam/H/buNShJkiQ1gcFIaiWLFsHxx+cNoFI5NCgNDBiUJEmSGsBgJLWy3l447ri8AYyP54BUDUsDA3mfJEmSjojBSGonPT1w7LF5gxyK9u7Na5MGB/M2NJRvK5Vya5UkSWojBiOpnfX0wOrVeas3MjIRlmoD0+BgbvogSZKkgwxGUqdasiRvk4WmsbHDw1J1Gx5ufq2SJEklMxhJ3WjRIli5Mm/1xscnD0xDQ3lzTZMkSepABiNJh+rpgaOOylu9lPKI0mTT81zXJEmS2pjBSNLsRcCyZXmrNoCoVV3XNNmIk+uaJElSCzMYSVo4061rqlQmn543MpJD09iY0/QkSVJpDEaSmqO3F1asyNtUKpUckEZHJ8LSdF9Xb73orSRJOkIGI0mto7c3b0uXzu28sbHZh6jq12NjBipJknSQwUhS+1u0KG/Lls3+nJSmDlTThayxscZ9H5IkqTQGI0ndKQIWL87b8uWzPy+liZBUqeR1UY3eHNmSJKnhGhaMIuItwLOBc4CRlNIxszgngMuB/wkcA3wTeEVK6bZG1SlJcxIx0WSiWRoVtqqBa7Lb2e6b72NTHS9JUkkaOWK0BPgscCPwP2Z5zp8Cfwy8FLgTeDtwfUT8UkppqCFVSlKr6+nJWzeYbchqFUdaz1Tnt+r+6s+gdptsfyscW7u/9hcCtb8gkKQaDQtGKaXLACLiZbM5vhgtugT4y5TSPxf7fg+4D/hN4FONqFOS1EIi8iY1w3ShaS73j+Tc+vuSStNKa4xOB9YBX6nuSCn1R8RNwFMwGEmSpIUUMdENs1XUh6bx8byesbqmcSG+brWRV6lFtFIwWlfc3le3/76axw4TEUuB2t6+Kxe4LkmSpOZoRlhLaeHDlsFLHWBOwSgi3gG8cYbDzkop/WT+Jc3ZpcBlTXw9SZKk9hUxcZmDRqgNXtUOntXbmfbVP+70QjXRXP9GvBv42AzH3DG/UthR3K4FttfsXwvcMs15VwBX1txfCdw9zxokSZJ0JGqD11wv2F2vGrKONGBVR7WkacwpGKWUdgI7G1TLneRw9EyKIBQRq4AnAX83TU3DwHD1frhoV5IkqTMs5OhW9cLeswlTtdMDZ7p16mDHaOR1jE4BjgNOAXoj4pzioZ+llPYVx/wEuDSl9IWUUoqIq4A/i4jbmGjXfS9wXaPqlCRJUheovbD3QqptknGkt1M95mhXUzSy+cLbyNcjqvpecXs+sKn4egOwuuaYvwaOBq4hX+D1G8CFXsNIkiRJLal6rblGrdmCiW6FCxGyprvtcpE6bPivmH7X39/fz6pVq8ouR5IkSWoP04WsCDj22LIrnLOBgQFWr14NsDqlNDDdsa3UrluSJElSWaqjX12qe79zSZIkSSoYjCRJkiR1PYORJEmSpK5nMJIkSZLU9QxGkiRJkrqewUiSJElS1+vYdt0DA9O2KZckSZLU4eaSCTrxAq8nAneXXYckSZKklnFSSume6Q7oxGAUwEOAvWXXMg8ryaHuJNqz/nbl+14O3/dy+L6Xw/e9HL7v5fB9L4fv+9RWAvemGYJPx02lK77hadNgq8qZDoC9KSXnAjaJ73s5fN/L4fteDt/3cvi+l8P3vRy+79Oa1fth8wVJkiRJXc9gJEmSJKnrGYxayzBweXGr5vF9L4fvezl838vh+14O3/dy+L6Xw/f9CHVc8wVJkiRJmitHjCRJkiR1PYORJEmSpK5nMJIkSZLU9QxGkiRJkrqewUiSJElS1zMYSZIkSep6BiNJkiRJXc9gJEmSJKnrGYwkSZIkdT2DkSRJkqSuZzCSJEmS1PUMRpIkSZK6nsFIktpIRKyNiM9FxO6ISBFxyQI851sjIi1AeTO9zmlFzX/S6NfS9I7kZx4RH4uIuxa4JEkqncFIktrLe4D/ClwB/C7w5TKLiYhXRsTLyqxhPiLit48kVEbEUUW4OG/BipIklWpR2QVIkubkGcA/p5T+puxCCq8EdgEfK7mOufpt4FHAVfM8/yjgsuLrTQtQjySpZI4YSVJ7OQHYU3YRkiR1GoORJDVQRJwaEe+PiK0RMVisDfpsRJw2ybHHRMR7IuKuiBiOiLsj4hMRsSYiXlasCQngVcVanSnXiNSu54mI10bEz4vX/3pEPGoWdS+KiD+PiNuLWu6KiL+KiKU1x9wFPBI4t1pPRGya5fsyZU0R8fvFcz1mkvPeHBGViDhxmudeGRFX1byP90fEv0fEY4vHNwHPBk6tqfuu4rElEfG2iPhORPRHxP6IuCEizq95/tOAncXdy2qe463F4+si4qPFz284IrZHxD/X/syL2v4lIs6LiJuL92FLdWpeRDy3uD9U1DLZe/GMorb9EbGneI2zJjnuVyLiP4vnuj0iLprmvfud4vUGI+KBiPhURJw81fE15724OG9vRAwUtb9mpvMkqZU4lU6SGusJwFOBTwF3A6cBrwA2RcQvpZQOAETECuAG4CzgI8B3gTXArwMnAZvJa4r+Afh34BOzfP3fA1YCfwssA14D/EdEnJ1Sum+a8z4MvBT4HPBu4EnApUV9zymOuQR4H7AP+F/Fvumec7Y1fa547CXA9+rOfQmwKaV0zzTP/wHg+cDVwI+A44FfKWr/blHravL7+trinH3F7SrgD4F/Aj5U1Pk/gOsj4okppVvIoegVwN8BXwCuLc69tbj9PDkwvg+4izzKdwFwSnG/6mHAPwIfBP4P8CfAFyPi5cBfAe8vjrsU+ExEbEgpjQNExLOALwF3AG8FlgOvBr4ZEY9NKd1VHHc28H+Lmt9K/n//cib5OUXEW4C3A58h//z7iufcHBGPSSntqT+nOO+C4v36KvDGYvdZwNOA9052jiS1pJSSm5ubm1uDNmD5JPueDCTgd2v2XV7se84kx0fN1wm4ehave1px7AHgxJr9Tyz2X1mz7635v4OD93+5OOZDdc/5rmL/+TX7fkAOKrN5L+ZS0z8C9wA9NfseUxz3shleZ89M7xHwL8Bdk+zvBZbU7TsG2AH8fc2+NUUtb53k2AT8yQyvf1dx3FNq9v2XmvfnlJr9f1TsP69m3/fI4ea4mn2PBirAx2v2fQEYrHu+s4Cxup/5qcW+N9fV+ShgtHY/eT3ZXTX3rwL6gd5m/b1yc3Nza8TmVDpJaqCU0mD164hYHBHHAz8jf3h/bM2hzwO+n1L6wiTPcSSttK9LNaMrKaVvAzcBvzbNOdXHrqzb/+7i9tlHUM9sa/oE8BDg/Jp9LyF/yP/8DM+/B3hSRDxkroWllCoppRGAiOiJiOPIoyw3c+jPayqDwAhwXkQcO8OxP0op3Vhz/6bi9j9SSr+YZP8ZRV3rgXOAj6WUHqip/VbyaOKvFcf1kjsYXlf7fCmlHwPX19XyXPL0+s9Enrq5JiLWkAPhbRz6c6i3BziaPComSW3LYCRJDRQRy4s1K9uAYXIHt53kkYXVNYc+lDz6stBum2TfT8mjN1M5FRgnB7iDUko7yB+CT21CTf8ObCeHISKiB/gtcke+vTM8/5+SRzq2RcS3I7fVPmO2xUXESyPiVmAI2E3+eT2bQ39ek0opDZOnk/0qcF9EbI6IP42IdZMc/ou6c/uLL7fVHVfdXw1a1fd/6yTP+WNgTUQcTZ4Kt5zJ3+/6cx9OXr92G/n7rd3OIk8HnMr7yT+/LxXrqj4SERdOc7wktSSDkSQ11vuAt5DXbbyQPF3qAvIH7lb/N7jhF32d8oVTqpCn0z0vIpaRRyweQl6LM9O5nyGPrrwauBd4A/DDiPjVmc6NiN8hTxW7nby26ELyz+s/mOXPK6V0FXAmeW3QEHndzo8naaBQmeIpptofs3n9eeoh/7yr32/9NmXDhpTS/eQRrF8H/j/yz+pLEfHxBtYrSQvO5guS1FjPJ6/5eH11R/FB/5i6424nj3IstIdPsu9MDm0CUO/n5A/KDyePQAAQEWvJdf+85tj5hKfZ1vQJ4PXAfyePwOzk8Clgk0opbSePZLw/Ik4gN114C7lhAUxd9/PJDQ2eWzuFMSIur3+JGV7/dvLUw3dHxMOBW4rv5XdmU/8Mqu//hkkeewSwK6W0PyKGyFP7Jnu/68+9nRy87kwp/XSuBRXTD79Ibh7RQ37vL4qIt6eUfjb92ZLUGlr9t5WS1O4qHP6b/leTF/nX+jzwyxHxnLr9RMSRjBT8Zm1r64h4IrnD3JemPoV/K24vqdv/uuL2X2v27efwkLcgNRVrZm4ld4l7HvCplNLYdE8cEb0RcciUt2JE415gac3u/Uw+Na46WnPwPY+IJwFPqTvuQHF7TN3rH1UE31q3A3vrXn/eitB3C/DSiDj4+pFbnv8Xip9fMep2Pfn9PqXmuLPIa49qXUv+3i+r//MW2fFT1VP/WMqd86od+hbke5akZnDESJIa61+A342IfnLr6KcAzyJPpav1LvJoxWcj4iPAd4DjyNOTXg58f56v/zPgGxHxd+QPqZcUr/3XU52QUvp+MQ3qj4oP3l8nd457KXkh/9dqDv8O8IqI+LPite5PKf3HAtb0CeBviq9nnEZHbq99d0R8jvye7SO/308gj9jU1v2iiLgS+E9gX0rpi+Sf13OBL0TEvwKnk9//HwErqienlAYj4kfFc/wUeIC8RmwR8NWI+Exxzhi5vflacsv2hfIGcpC8MSL+nol23f3kLoNVl5Gnx90QEe8v6ns18ENyF7vq93N78TO8AjgtIq4jh7nTi/qvYeLnUO/DRZOK/yC3pD+1eI1bqBlxlKSWV3ZbPDc3N7dO3sgjCh8hTwPbC3yZPI3pLnJXsdpjjyOvSbqb3KhhG3m9y/E1x8y1XfefkEd6fkFe77IZeHTdsW+lpnVzsW8R8BfkaWUjxfl/BSytO24tOUwMFK+3aSFqqjlnHTlcbJ3l+72EHLBuKWraV3z9irrjjgY+CTxY1HRXsT/Ia4PuKmr7Lrnxwseoa+9NDrk3Fz+rVLyPx5Ovn/Tj4rX3AP8PeEHduXcB/zJJ/Yf9fGvft7r9zwS+QR696iev7zlrkufcWFPn7eT1Qof9zItjn0u+nta+Yvtx8f2cWXPMIe8FeTTvenL78GHyVL8PAOvK/vvn5ubmNpctUiptba0kqUEi4jTgTuANKaWpftPf8oqW0duBt6WU3l52PZKkzuUaI0lSK3sZeT3WP5RchySpw7nGSJLUciLiGcAvkTvJXZdSuqvciiRJnc5gJElqRX8BPBX4JnkhvyRJDeUaI0mSJEldzzVGkiRJkrqewUiSJElS1+u4NUbFFbsfQr5eiCRJkqTuthK4N82whqjjghE5FN1ddhGSJEmSWsZJwD3THdCJwWgvwLZt21i1alXZtUiSJEkqycDAACeffDLMYjZZJwYjAFatWmUwkiRJkjQrNl+QJEmS1PUMRpIkSZK6nsFIkiRJUtczGEmSJEnqegYjSZIkSV3PYCRJkiSp63Vsu25JkiRJc1MZT4ynRGU8kRJUiq+XL+5lyaLOHlNpaDCKiI3AG4DHAeuB56SUrpvhnPOAK4FHAtuAv0wpfayRdUqSJEmtLBUBZTzBeJoIL+MJxqthJhVhprg/Pp6DTf46H3vwsTRxv/a5p/LQE47mhJXLmvcNl6DRI0ZHA98HPgJcO9PBEXE68K/AB4CXAM8EPhwR21NK1zey0LLduWs/n7l5G3c/OMhJxy7nhY8/mdPXHF12WZIkSZrBnILJwa/rgklx/GGhp2b0Ro3V0GCUUvoS8CWAiJjNKS8H7kwpvb64/+OI+BXgtUDHBqPP3LyNN33+ViKClBIRwQe/fjvvfN6jecHjTy67PEmSpI5QGU8TW0qH3h9Ph4SWSl14Ga8JLjkAJSrFfUNLZ2i1NUZPAb5St+964Krml9Icd+7az5s+f2seuqz+rSpu3/j5W3nCacdxmiNHkiSpS00XZsZTYmw8h5mx2n2ViXAzVrdmRppKqwWjdcB9dfvuA1ZFxPKU0mD9CRGxFFhas2tlA+tbcJ+5eVseTZvkb2pE8Ombt/HGCx9RQmWSJEnzUx9cKjXhpTbM1AeXyiT7DDNqllYLRvNxKXBZ2UXM190PDpKm+BufUuLuBw/LgpIkSUdkvDptbJJ1MPVrWyY9pnZaWc1ifqeWqZ21WjDaAayt27cWGJhstKhwBbmLXdVK4O4G1NYQJx27fNoRo5OOXV5CVZIkqSwpHdp5zOAiNUerBaMbgV+r23dBsX9SKaVhYLh6f5ZNHlrGCx9/Mh/8+u2TPpZS4kU2X5AkqWFSERLGa8IITNxPM90y0S451d4vFufn554II4fcMnFcfk5Di1SmRl/HaAXwsJpdp0fEOcADKaVfRMQVwIkppd8rHv8AcHFE/DW5xfczgBcCz25knWU6fc3RvPN5j+aN1QYMQG8EicQ7n/doGy9I6nqp7oNiO39mnGzq9FTfz1QfjtMUZ8z1w/Rcnj9/4K9+nY9ICUgTx1ePqT4+sS8Vx9X0GCIdcnxtPdVzJo6te83qM9Q85yE11J8/i7AiSdD4EaPHA1+ruV+d8vZx4GXki76eUn0wpXRnRDwbeA/wGvKUuD/s9GsYveDxJ/OoE1fxq+/9BgC//yun8TtPOtVQJOmg+qk1h1wnI03++MHfgI/X3q/9gFn3gfTgax36QXOqfUx3PDUfeGuOqz462WtWv8/Dj5ckqfEafR2jTcCUc9tSSi+b4pzHNKyoFnXq8RMh6HUXnMlRS1ptlqMkyB/cqx2WqmGjkhKp5qJ8MwaVmqkztdfFqA0v9c9hSJAkqbH89C2pa41Vxhkr2sIe/LqSGK2MF+1iJ9tnSJEkqRMZjCS1vYMhppJDTDXQVMPMaKXm8fFEZXyc0YoBR5IkTTAY6YjduWs/n7l5G3c/OMhJxy7nhY8/mdNdH6V5GB9PjI5PhJlKMZIzOp6oVPJjY4eEnPHi6uZlVy5JktqdwUhH5DM3b+NNn7+ViCClRETwwa/fzjuf92heYKtx1RirjDM0Ns7QaKXYxhkZG2esGL2pXu1ckiSpDAYjzdudu/bzpmqb8drWVcAbP38rTzjtODvrdZnhsRx4hscqDI+OHwxAQ2MVxiqGHkmS1LoMRpq3z9y8LV9Qd5KFGhHBp2/exhsvfEQJlalRUkoMj9UEntFKzf2KU9okSVLbMhhp3u5+cHDSixVC/gB994ODTa5IC6Eyng6O/NROexsaqzAyNm7DAkmS1JEMRpq3k45dPu2I0UnHLi+hKs3GaGVi1Kc2BA2PVRgZM/lIkqTuYzDSvL3w8Sfzwa/fPuljKSVeZPOFUh1c71MfgFzvI0mSdBiDkebt9DVH887nPZo3VhswAL0RJBLvfN6jbbzQBKOVcfYPj01Me6sJQ673kSRJmj2DkY7ICx5/Mo86cRW/+t5vAPD7v3Iav/OkUw1FDTIyNs7eoVEGhsYYGBzlwEil7JIkSZI6gsFIR+zU4ydC0OsuOJOjlvjHaqEMj1XYW4SggaExBg1CkiRJDeEnWKmFDI1WGBgaPRiGhkbHyy5JkiSpKxiMpBINjVYOjgYNDI0ybBCSJEkqhcFIaqLBkeqI0Cj9g2OMjBmEJEmSWoHBSF3tzl37+czN27j7wUFOOnY5L3z8yZy+gI0jDoyMMTA4djAMeY0gSZKk1mQwUtf6zM3beNPnbyUiSCkREXzw67fzzuc9mhfM4xpMKSUOFCNCA4Nj7B0aZdTrBUmSJLUFg5G60p279vOm6vWXUhFeits3fv5WnnDacTO2HE8psX+kukYoN0zwwqmSJEntyWCkrvSZm7cREROhqEZE8Ombt/HGCx9xyP6UEvuGxw5eQ2jv0BgVr6IqSZLUEQxG6kp3PzhImiQUQQ5Adz84yPh4Yu9wnhI3MDjGvmGDkCRJUqcyGKkrnXTs8ilHjACWLerhP+96AHOQJElSd+gpuwCpDC98/MlTjxgBT33oGkORJElSFzEYqevsGx6jN4JXnf8wIib29wREwEUbz2Dd6mXlFShJkqSmcyqdOl5KiYHBMR44MMID+0cOXlT1qQ9dw4nHLOdN124B4MJHreOCs9YZiiRJkrqQwUgdaXw8sWdwlAf2j7DnwMiU1xNau2oiBL3gcSezbHFvs0qUJElSCzEYqWOMVcZ54MAID+4fpX9w1A5ykiRJmjWDkdra8FiFB/fnkaGBodGpmsxJkiRJ0zIYqe0MjlSKkaER9g6NlV2OJEmSOoDBSG1h79BoHhk6MMLgSKXsclrC9v5BNm3dyc59w/StWMp5G/pYv3p52WVJkiS1JYORWtJUneSUbdp6P9fccAdBvu5SAF+89V4u2ngG5555QsnVSZIktR+DkVpGZTzRPzjKA/uHefDAKGNTdJLrdtv7B7nmhjtIKYcimLj94OY72LB2lS3HJUmS5shgpFKNVsZ5sOgkt+fACDaSm9mmrTsPjhTVC+BrW+/nt554SpOrkiRJam8GIzWdneSOzM59w5OGIshhaee+4WaWI0mS1BEMRmqKwZEKu/cP8+D+UfYN20nuSPStWDrtiFHfiqVNrkiSJKn99TTjRSLiVRFxV0QMRcRNEfHEaY59WUSkum2oGXVqYe0dGuUXuw9wy7Y93LJtD9seGDQULYDzNvRNO2J0/gabL0iSJM1Vw0eMIuJFwJXAy4GbgEuA6yNiQ0rp/ilOGwA21Nx3slWb6D8wyn0Dw3aSa6D1q5dz0cYz+ODmOw5OQ+yJ/Jfkoo1n2HhBkiRpHpoxle51wIdSSh8FiIiXA88G/gB4xxTnpJTSjibUpiNUGU88sG/k4P2f7NjLssW9JVbUHc498wROO/5o3nTtFgAufNQ6LjhrnaFIkiRpnhoajCJiCfA44IrqvpTSeER8BXjKNKeuiIifk6f6fRd4c0rph1O8xlKgdlHFyiMuXNOq7yR3wAuulmLtqokQ9ILHnWwglSRJOgKNHjFaA/QC99Xtvw94xBTnbCWPJt0KrAb+BPhWRDwypXT3JMdfCly2MOVqKnaSkyRJUidrua50KaUbgRur9yPiW8CPgYuAP5/klCvIa5iqVgKTBSjNkZ3kJEmS1C0aHYx2ARVgbd3+tcCs1hCllEYj4nvAw6Z4fBg4eOGWiJhfpQJyJ7kH94/ywIERBp0iJ0mSpC7R0GCUUhqJiO8AzwSuA4iInuL+1bN5jojoBc4G/q1BZXa1lBIDg2N5ZOjAqJ3k1BTb+wfZtHUnO/cN07diKedt6GP96uVllyVJkrpYM6bSXQl8PCJuBr5Nbtd9NFDtUvcJ4J6U0qXF/b8A/h/wM+AY4A3AqcCHm1BrV6iMJ/YcGMkNFA6MMlZxwZCaZ9PW+7nmhjsOXqQ2gC/eei8XbTyDc8/0GkySJKkcDQ9GKaVPR0Qf8DZgHXALcGFKqdqQ4RSgdpjiWOBDxbEPAt8BnppS+lGja+10I2Pj3LlrP3sOjDBuFlIJtvcPcs0N+fpL1T+C1dsPbr6DDWtX2XJckiSVoinNF1JKVzPF1LmU0nl1918LvLYJZXWdkco4D+wfmflAqUE2bd15cKSoXgBf23o/v/XEU5pclSRJUr5OkCQ1xc59w5OGIshhaee+4SkelSRJaiyDkaSm6VuxlKn6RkbxuCRJUhkMRpKa5rwNfdOOGJ2/weYLkiSpHAYjSU2zfvVyLtp4BrWXG+sJiICLNp5h4wVJklSapjRfkKSqc888gdOOP5o3XbsFgAsftY4LzlpnKJIkSaUyGElqurWrJkLQCx53MssW95ZYjSRJklPpJEmSJMlgJEmSJElOpZOkOdjeP8imrTvZuW+YvhVLOW9DH+tXLy+7LEmSdIQMRpI0S5u23s81N9xBkNuLB/DFW+/loo1ncO6ZthqXJKmdOZVOkmZhe/8g19xwBynBeOKQ2w9uvoMd/UNllyhJko6AwUiSZmHT1p3EFI8F8LWt9zezHEmStMAMRpI0Czv3DZOmeCwVj0uSpPZlMJKkWehbsXTaEaO+FUubWY4kSVpgBiNJmoXzNvRNO2J0/gabL0iS1M4MRpI0C+tXL+eijWcQNcNGPQERcNHGM1i3ell5xUmSpCNmu25JmqVzzzyB044/mjdduwWACx+1jgvOWmcokiSpAxiMJGkO1q6aCEEveNzJLFvcW2I1kiRpoTiVTpIkSVLXMxhJkiRJ6npOpZOkLrG9f5BNW3eyc98wfSuWct6GPtavXl52WZIktQSDkSR1gU1b7+eaG+4gyO3FA/jirfdy0cYzOPdMW41LkuRUOknqcNv7B7nmhjtICcYTh9x+cPMd7OgfKrtESZJKZzCSpA63aetOYorHAvja1vubWY4kSS3JYCRJHW7nvmHSFI+l4nFJkrqdwUiSOlzfiqXTjhj1rVjazHIkSWpJBiNJ6nDnbeibdsTo/A02X5AkyWAkSR1u/erlXLTxDKJm2KgnIAIu2ngG61YvK684SZJahO26JakLnHvmCZx2/NG86dotAFz4qHVccNY6Q5EkSQWDkSR1ibWrJkLQCx53MssW95ZYjSRJrcWpdJIkSZK6niNGkqSWt71/kE1bd7Jz3zB9K5Zy3oY+1q9eXnZZkqQOYjCSJLW0TVvv55ob7iDIXfQC+OKt93LRxjM490w76kmSFkZTptJFxKsi4q6IGIqImyLiiTMc/4KI+Elx/JaI+LVm1ClJai3b+we55oY7SAnGE4fcfnDzHezoHyq7RElSh2h4MIqIFwFXApcDjwW+D1wfEZP+mi8ingr8E/D3wGOA64DrIuJRja5VktRaNm3dOe3Fab+29f5mliNJ6mDNmEr3OuBDKaWPAkTEy4FnA38AvGOS418DfDml9K7i/p9HxAXAxcDLZ/uiB0bGWDQydkSFN9OBmloPNKjuAyNjDI1WFvx5h2uec7gBz99I1l4Oay9HO9Z+38DQtBenvW9gqCH/rkmSDjU4UmnYZ9RGmkvNkdJU/+UcuYhYAhwAnp9Suq5m/8eBY1JKvzHJOb8ArkwpXVWz73LgN1NKvzzJ8UuBpTW7VgJ3n3zJZ+hZetRCfSuSJEmS2sz48AG2XfVCgNUppYHpjm30VLo1QC9wX93++4B1U5yzbo7HXwr012x3z6tSSZIkSV2rE7rSXUFew1S1Erj72295JqtWrSqppNa0b3iMH94zbVCWpJbzjdt28ZFv3XlIV7oE/MFTT+dXHr6m3OJmYXi0wss/+V0APvCSx7K0jS6sa+3lsPZyWPv0zug7mr6VS2c+sMUMDAyw/qrZHdvoYLQLqABr6/avBXZMcc6OuRyfUhoGhqv3I/Iy3aOWLOKoJZ2Q+xbOeMIr3UtqO8/6pbU86sTVfG3r/QevY3T+hhNYt3pZ2aXN2dLFvW3777C1l8Pay2Hth1u+pLctP1uPzaHmhn53KaWRiPgO8Exydzkioqe4f/UUp91YPH5Vzb4Liv2SpC60bvUyfuuJp5RdhiSpgzXjOkZXAv8zIl4aEWcBfwccDVS71H0iIq6oOf69wIUR8fqIeEREvBV4PFMHKUmSWtaOgYlrLX32O9vY3j9YYjWSpKk0fDwspfTpiOgD3kZuoHALcGFKqdpg4RRgvOb4b0XEbwN/CfwVcBu5I90PGl2rJEkLadPW+7nmhjsO3v/yD3bwpR/s4KKNZ3DumZNezk+SVJKmTBRMKV3NFCM+KaXzJtn3WeCzDS5LkqSG2d4/yDU33EHtVTHGi68/uPkONqxd1ZbrpCSpUzVjKp0kSV1n09adxBSPBfC1rfc3sxxJ0gwMRpIkNcDOfcNMdQn1VDwuSWodBiNJkhqgb8XSaUeM+la03/VAJKmTGYwkSWqA8zb0TTtidP4Gmy9IUisxGEmS1ADrVy/noo1nEAE9wSG3F208w8YLktRi2u/ytZIktYlzzzyBDWtX8bWt97Nz3zB9K5Zy/oYTDEWS1IIMRpIkNdC61cv4rSeeUnYZ81J/cdpnnbWW9auXl1iRJDWOU+kkSdJhNm29nzd/YcvB+1/+wQ5e/9nv8/Wf2mZcUmcyGEmSpENMdXHalPLFaXf0D019siS1KYORJEk6hBenldSNDEaSJOkQXpxWUjcyGEmSpEN4cVpJ3chgJEmSDuHFaSV1I4ORJEk6RCdcnLa+1fj2/sESq5HUDryOkSRJOkw7X5x209b7ueaGOw7e//IPdvClH+zgoo1ncO6ZjnZJmpzBSJIkTaodL047VatxyK3GN6xd1RbhTlLzOZVOkiR1DFuNS5ovg5EkSeoYthqXNF8GI0mS1DFsNS5pvgxGkiSpY3RCq3E76knlMBhJkqSO0e6txjdtvZ83f2HLwftf/sEOXv/Z7/P1n7o2Smo0u9JJkqSO0q6txu2oJ5XLYCRJkjpOO7Yar3bUm2wqYLWjXrt9T1I7cSqdJElSC+iEjnquj1I7MxhJkiS1gHbvqOf6KLU7g5EkSVILaOeOelOtj0opr4/a0T809clSizAYSZIktYB27qhXXR81mer6qFbnNEDZfEGSJKlFtGtHvXZfH7Vp6/1cc8MdB+9/+Qc7+NIPdnDRxjM498zWHamrqg91zzprLetXLy+xovZkMJIkSWoh7dhRr7o+aqqOeq28Pqrd26S3e6hrJU6lkyRJ0hFp5/VR7TwN0LVdC8tgJEmSpCPSzuuj2nkaYDuHulbkVDpJkiQdsXZdH9XO0wDbOdS1IoNRFzlqcS9nrl3BgwdGePDAKGOVqf4qSZIkzV07ro86b0MfX7z13kkfa/VpgO0c6lpRQ6fSRcRxEfHJiBiIiD0R8fcRsWKGczZFRKrbPtDIOrtFT09w/IqlPOyElTzulGM5a/1K1q5aypJFUw3CSpIkdbZ2ngbYzmu7WlGjR4w+CawHLgAWAx8FrgF+e4bzPgT8Rc39Aw2prov19ATHHLWEY45aQkqJvcNjPLh/hAf2jzA0Ol52eZIkSU3TrtMAq6Hug5vvODhyVL1t9VDXiiKlxkynioizgB8BT0gp3VzsuxD4N+CklNKkY5YRsQm4JaV0yTxfdxXQ39/fz6pVq+bzFF3vwMgYDxQhaf9wpexyJEmSNI0d/UMND3UPPeFoTljZfkFrYGCA1atXA6xOKQ1Md2wjR4yeAuyphqLCV4Bx4EnAF6Y59yUR8TvADuCLwNtTSo4aNclRSxZx1JJFnHTsUQyNVnjwQA5Je4fGaFCOliRJ0jy149quVtTIYLQOOKRHYEppLCIeKB6byj8CPwfuBR4NvBPYADx3soMjYilQu7Js5RHUrDrLFveyfvVy1q9ezmhlPE+3OzBC/4HRgxc/kyRJktrdnINRRLwDeOMMh501v3IgpXRNzd0tEbEd+GpEPDSldPskp1wKXDbf19PsLe7t4YRVyzhh1TIq44k9B0bscCdJkqSOMJ8Ro3cDH5vhmDvI0+AOaYUREYuA44rHZuum4vZhwGTB6Argypr7K4G75/D8mofeosPd8SuWklKif3CUB/bnoDQyZkiSJElSe5lzMEop7QR2znRcRNwIHBMRj0spfafY/Qxyi/Cbpj7zMOcUt9unqGcYOHj1qghbTzdbxESHO4C9Q6MHmzfY4U6SJEntoGFrjFJKP46ILwMfioiXk9t1Xw18qtqRLiJOBL4K/F5K6dsR8VByK+9/A3aT1xi9B9icUrq1UbVqYa1ctpiVyxZz6vFHH+xw9+D+UfYNj5VdmiRJkjSpRl/H6CXkMPRVcje6zwN/XPP4YnJjhaOK+yPAs4BLgKOBbcU5f9ngOtUgEx3usMOdJEmSWlbDrmNUFq9j1B5GK+MHQ5Id7iRJklqb1zGSGmRxbw8nrFzGCSvtcCdJkqTyGYxUOjvcSZIkqWwGI7WUyTrcPbh/lN37h+1wJ0mSpIYxGKmlVTvcnXL8UXa4kyRJUsMYjNQ27HAnSZKkRjEYqS0tW9zL+tXLWb96uR3uJEmSdMQMRmp7driTJEnSkTIYqaPY4U6SJEnzYTBSx7LDnSRJkmbLYKSuYYc7SZIkTcVgpK5U3+Guf3CUgcFRBoZGnXInSZLUhQxG6nrLFveybHEva1ctA2BwpMLeoRyS+gfHGBlz2p0kSVKnMxhJdZYv6WX5kl5OKILS0GiFgaFRBgbHGBgaZdj1SZIkSR3HYCTNoDqidMLKfH94rHIwJA0MjtrIQZIkqQMYjKQ5Wrqol76VvfStXArAyNj4wZC0d2iMAyOVkiuUJEnSXBmMpCO0ZFEPa1YsZc2KHJRGK+NFI4cx9g6Nsn/YoCRJktTqDEbSAlvc23PwIrMAY5XxgyFpYHCM/SNjJBvfSZIktRSDkdRgi3p7OO7oJRx3dL7Q7FhlnL1DY+wdyuuU9g0blCRJkspmMJKabFFvD8cevYRji6BUGU/sK0JS/+Ao+4fHGDcoSZIkNZXBSCpZb0+w+qjFrD5qMScD4+OJvcNjBy84u2/IoCRJktRoBiOpxfT0BKuXL2b18sVADkr7RoqgNDjGvuExKiYlSZKkBWUwklpcT0+watliVi1bDMdCSol9w2MHGzrsGxpjtGJQkiRJOhIGI6nNRAQrly1m5bLFwHIgN3QYGhtnaLTC0GiF4YNfjzMy5gVoJUmSZmIwkjrAot4eVvT2sGLp4X+lx8cTQ2MVhkfHGRrLYak2QNkRT5IkyWAkdbyenuCoJYs4asnhj6WUGB4brwlNOTgNFwHKtUySJKlbGIykLhYRLFvcy7LFvaxm8WGPj4xNBKbhmsA0NFpxXZMkSeooBiNJU1qyqIcli3py44c61XVNw6OVQ9Y3ua5JkiS1I4ORpHmZaV3TwQYQNaNMw0WQcoaeJElqNQYjSQuupydYvqSX5Ut6D3usdl3TcKXCWCXlbXycsfHEaCWvbRqtJCrjyXVOkiSpKQxGkpqqdl0Tk6xrqjc+nhgdnwhLY9XgNJ6/HhuvCVY1t+YpSZI0FwYjSS2tpydY2nP4yNNMKjWjT2OV2nBVE7IOhqmJkGX7ckmSupPBSFJH6u0JeucRqA6OQtWNSB0MWePjB6f5pQTjKVFJiZTyKNX4uKNVkiS1I4ORJNVY1NvDornnqcPkgFSEpTQRoqr7Us1jhz6ez53s+EpKjI/P/HyOekmSNHcNC0YR8Rbg2cA5wEhK6ZhZnBPA5cD/BI4Bvgm8IqV0W6PqlKRG6OkJeohSXnva0JUSqeimnkjF8ZBqzs2PUfNYonbnZI+lwx47PKBNHJMmP77uOCap79DHW0difkVN9b1MuX+K15nrezKf56//WeXnOPRnmdJUP0tJan2NHDFaAnwWuBH4H7M850+BPwZeCtwJvB24PiJ+KaU01JAqJanDRAS9Ab0lBTOpVjUkJw4N3tVgXRueavdNGsTqwnXi8CCWyA9U99X+cqD+drwmyI3XjMhWz6sdhXVUVup8DQtGKaXLACLiZbM5vhgtugT4y5TSPxf7fg+4D/hN4FONqFOSJDVORBAHM3rnhPXqdNlqqJoydE0VwoppszB9CKuuZ6wc/DrfupZRWnittMbodGAd8JXqjpRSf0TcBDwFg5EkSWoRZU6XhYnwVDm4nrEITOMTTWGqzWAmvk41x3PwWnEGLylrpWC0rri9r27/fTWPHSYilgJLa3atXOC6JEmSWsrBKbM9jQlnswpexfTDw4PXROgaT7nLZ6XmcaciqlXNKRhFxDuAN85w2FkppZ/Mv6Q5uxS4rImvJ0mS1NEaGbyqlz6ohqpKJQetyfZVA9bBrW6ftJDmOmL0buBjMxxzx/xKYUdxuxbYXrN/LXDLNOddAVxZc38lcPc8a5AkSVIDzfc6c5MZq4wfvJRBbbCqfn1I2BofpzLFPkOWYI7BKKW0E9jZoFruJIejZ1IEoYhYBTwJ+LtpahoGhqv3IzpnYackSZKmtqi3p+bD7PzDVnVtVaV2emDtOq3q+quaxhmTrdmq7W5Yvfh3Zdxuhu2ikdcxOgU4DjgF6I2Ic4qHfpZS2lcc8xPg0pTSF1JKKSKuAv4sIm5jol33vcB1japTkiRJ3S0iWNQbDV18fzBMFSFpsiA1WcOMw4LZOIeu+6p2PrRpxhFr5M//beTrEVV9r7g9H9hUfL0BWF1zzF8DRwPXkC/w+g3gQq9hJEmSpHZW7WTYyA/fUzbNSLMLZrUjXPUdC3u6YFZWpA4b1yum3/X39/ezatWqssuRJEmSVJKBgQFWr14NsDqlNDDdsT3NKUmSJEmSWpfBSJIkSVLXMxhJkiRJ6noGI0mSJEldz2AkSZIkqesZjCRJkiR1PYORJEmSpK7XyGtMlWpgYNo25ZIkSZI63FwyQSde4PVE4O6y65AkSZLUMk5KKd0z3QGdGIwCeAiwt+xa5mElOdSdRHvW365838vh+14O3/dy+L6Xw/e9HL7v5fB9n9pK4N40Q/DpuKl0xTc8bRpsVTnTAbA3peRcwCbxfS+H73s5fN/L4fteDt/3cvi+l8P3fVqzej9sviBJkiSp6xmMJEmSJHU9g1FrGQYuL27VPL7v5fB9L4fvezl838vh+14O3/dy+L4foY5rviBJkiRJc+WIkSRJkqSuZzCSJEmS1PUMRpIkSZK6nsFIkiRJUtczGLWQiHhVRNwVEUMRcVNEPLHsmjpZRFwaEf8ZEXsj4v6IuC4iNpRdV7eJiDdFRIqIq8qupdNFxIkR8X8iYndEDEbEloh4fNl1dbKI6I2It0fEncV7fntE/HnUXIlRRy4iNkbEFyPi3uLfk9+sezwi4m0Rsb34OXwlIh5eUrkdY7r3PSIWR8Q7i39n9hfHfCIiHlJiyR1hpj/vdcd+oDjmkuZV2L4MRi0iIl4EXElus/hY4PvA9RFxQqmFdbZzgb8FngxcACwG/m9EHF1qVV0kIp4AXATcWnYtnS4ijgW+CYwCvwr8EvB64MEy6+oCbwReAVwMnFXc/1Pg1WUW1YGOJv+/+aopHv9T4I+BlwNPAvaT/49d1pzyOtZ07/tR5M8zby9unwtsAP6/plXXuWb68w5ARDyH/Bnn3mYU1Qls190iIuIm4D9TShcX93uAbcD7UkrvKLW4LhERfcD9wLkppc1l19PpImIF8F3glcCfAbeklC4ptagOFhHvAJ6WUnp62bV0k4j4F+C+lNL/qNn3eWAwpfQ75VXWuSIiAc9JKV1X3A/yB8N3p5T+pti3GrgPeFlK6VNl1dpJ6t/3KY55AvBt4NSU0i+aVVsnm+p9j4gTgZuA/wr8K3BVSumqphfYZhwxagERsQR4HPCV6r6U0nhx/yll1dWFVhe3D5RaRff4W+BfU0pfmfFILYRfB26OiM8WU0e/FxH/s+yiusC3gGdGxJkAEfHLwK8AXyq1qu5yOrCOQ/+P7Sd/aPT/2OZaDSRgT8l1dLTil+v/ALwrpfTDsutpJ4vKLkAArAF6yb+9qnUf8Ijml9N9in9ErgK+mVL6QcnldLyIeDF5asUTyq6li5xBntJ1JfBX5Pf+f0fESErp46VW1tneAawCfhIRFfK/9W9JKX2y3LK6yrridrL/Y9ehpiimLb4T+KeU0kDZ9XS4NwJjwP8uu5B2YzCSsr8FHkX+Ta4aKCJOBt4LXJBSGiq7ni7SA9ycUnpzcf97EfEo8poLg1HjvBB4CfDbwA+Bc4CrIuJeA6m6RUQsBj4DBPkXNGqQiHgc8Brgscn1MnPmVLrWsAuoAGvr9q8FdjS/nO4SEVcD/w04P6V0d9n1dIHHAScA342IsYgYIzfC+OPifm+55XWs7cCP6vb9GDilhFq6ybuAd6SUPpVS2pJS+gfgPcClJdfVTar/j/p/bAlqQtGp5F+IOVrUWE8n/x/7i5r/Y08F3h0Rd5VaWRswGLWAlNII8B3gmdV9xdSuZwI3llVXpyvat14NPAd4RkrpzrJr6hJfBc4m/+a8ut0MfBI4J6VUKauwDvdNckeoWmcCPy+hlm5yFDBet6+C//82053kAFT7f+wqcnc6/49toJpQ9HDgWSml3SWX1A3+AXg0h/4fey/5lzT/tayi2oVT6VrHlcDHI+JmcseWS8jtGD9aZlEd7m/J01t+A9gbEdW55v0ppcHyyupsKaW9wCHruCJiP7Db9V0N9R7gWxHxZvIHlScCf1RsapwvAm+JiF+Qp9I9Bngd8JFSq+owRZfLh9XsOj0izgEeSCn9orhO2p9FxG3koPR28ofF65pcakeZ7n0nj1J/jrye9L8BvTX/zz5Q/FJY8zDTn3dgd93xo8COlNLW5lXZnmzX3UIi4mLgDeTFoLcAf5xSuqnUojpY0eJyMr+fUvpYM2vpdhGxCdt1N1xE/DfgCvJvb+8ErkwpfajcqjpbRKwkfwh/Dnl6y73APwFv84PhwomI84CvTfLQx1NKLytadl9O/kXAMcA3gFemlH7arBo70XTvO/BW8r8zkzk/pbSpIUV1gZn+vE9y/F3YrntWDEaSJEmSup5znCVJkiR1PYORJEmSpK5nMJIkSZLU9QxGkiRJkrqewUiSJElS1zMYSZIkSep6BiNJkiRJXc9gJEmSJKnrGYwkSZIkdT2DkSRJkqSuZzCSJEmS1PUMRpIkSZK63v8POu/5H2my5AsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 自己画一个\n",
    "# 自己计算\n",
    "acf_value, acf_interval, _, _ = acf(data.value,nlags=14,qstat=True,alpha=0.05, fft=False)\n",
    "\n",
    "xlabel = np.arange(start=0, stop=acf_value.shape[0], dtype='float')\n",
    "\n",
    "fig, ax = plt.subplots(nrows=2, figsize=(10,6), sharex=True, dpi=100)\n",
    "ax[0].hlines(y=0, xmin=np.min(xlabel)-2, xmax=np.max(xlabel)+2, colors='red', linewidth=0.5)\n",
    "ax[0].scatter(x=xlabel, y=acf_value, c='red')\n",
    "ax[0].vlines(x = xlabel, ymin=0, ymax=acf_value, colors='red', linewidth=0.5)\n",
    "xlabel[1] -= 0.5\n",
    "xlabel[-1] += 0.5\n",
    "ax[0].fill_between(x=xlabel[1:], y1=acf_interval[1:,0] - acf_value[1:], y2=acf_interval[1:, 1]-acf_value[1:], alpha=0.25, linewidth=0, color='red')\n",
    "ax[0].set_title(\"acf plot by myself\")\n",
    "\n",
    "\n",
    "# 使用别人写好的函数\n",
    "\n",
    "from statsmodels.graphics.tsaplots import plot_acf\n",
    "plot_acf(data['value'], ax=ax[1])\n",
    "ax[1].set_title(\"acf plot by statsmodels\")\n",
    "ax[1].set_xlim(-1, np.max(xlabel)+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# 偏自相关系数\n",
    "\n",
    "偏自相关系数可以被看成一种条件相关，两个变量之间有相关性是基于一系列的其他的条件变量。\n",
    "想象一个这样的问题：\n",
    "\n",
    "我们有一个Y变量，这个Y变量和$X_1$、$X_2$相关，还有一个变量$X_3$。那我要知道$Y$和$X_3$是否相关应该怎么计算？\n",
    "\n",
    "肯定是应该要剔除掉$X_1$、$X_2$对$Y$、对$X_3$的影响。那么计算$Y$和$X_3$的相关性就是为：\n",
    "\n",
    "$\\gamma(X_3, Y) = \\frac{Cov(X_3, Y | X_1,X_2)}{\\sqrt{Var(X_3|X_1,X_2)Var(Y|X_1, X_2)}}$\n",
    "\n",
    "那么对于一个时间序列来说，$X_t$ 和$X_{t-h}$之间的偏自相关系数：就是基于 $X_{t-h+1}$、$X_{t-h+2}$、$X_{t-h+3}$、…… $X_{t-1}$对$X_t$ 和$X_{t-h}$的影响，计算$X_t$ 和$X_{t-h}$之间相关系数。\n",
    "\n",
    "## 偏自相关计算公式：\n",
    "1. 当h=1的时候，我们考虑的就是$X_t$ 和$X_{t-1}$之间的偏自相关系数,这个时候中间是没有别的条件序列的，所以偏自相关系数就是等于自相关系数。\n",
    "2. 当h=2的时候，我们考虑的就是$X_t$ 和$X_{t-2}$之间的偏自相关系数，这个时候中间的条件序列就是$X_{t-1}$，所以偏自相关系数就是为:\n",
    "\n",
    "$\\gamma(X_{t-2}, X_t) = \\frac{Cov(X_{t-2}, X_t| X_{t-1})}{\\sqrt{Var(X_{t-2}|X_{t-1})Var(X_{t}|X_{t-1})}}$.\n",
    "\n",
    "3. 当h=3的时候，我们考虑的就是$X_t$ 和$X_{t-3}$之间的偏自相关系数，这个时候中间的条件序列就是$X_{t-1}$、$X_{t-2}$，所以偏自相关系数就是为:\n",
    "\n",
    "$\\gamma(X_{t-3}, X_t) = \\frac{Cov(X_{t-3}， X_t | X_{t-1}，X_{t-2} )}{\\sqrt{Var(X_{t-3} | X_{t-1}，X_{t-2} ) Var(X_{t} | X_{t-1}，X_{t-2} )}}$.\n",
    "\n",
    "4. 依次类推。\n",
    "\n",
    "## 偏自相关计算：\n",
    "\n",
    "1. 使用statsmodels的pacf函数，这个函数可以选择不同的计算方法，有使用ols方法的，也有使用Yule-Walker方法的。\n",
    "2. 除了方法1调用包，我们还可以使用公式方法，分别使用ols方法、使用Yule-Walker方法来计算pacf的系数。实际上算出来的结果和statsmodels的pacf函数算出来的结果是一模一样的。\n",
    "\n",
    "在下面的代码块中，我们使用statsmodels包里面的pacf、pacf_ols、pacf_yw函数计算偏自相关系数，另外我们自己通过ols理论和Yule-Walker理论计算偏自相关系数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pacf</th>\n",
       "      <th>pacf_ols</th>\n",
       "      <th>pacf_yw</th>\n",
       "      <th>pacf_ols_bymyself</th>\n",
       "      <th>pacf_yw_bymyself</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.997827</td>\n",
       "      <td>0.997827</td>\n",
       "      <td>0.904116</td>\n",
       "      <td>0.997827</td>\n",
       "      <td>0.904116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.700485</td>\n",
       "      <td>-0.700485</td>\n",
       "      <td>-0.122213</td>\n",
       "      <td>-0.700485</td>\n",
       "      <td>-0.122213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.051486</td>\n",
       "      <td>-0.051486</td>\n",
       "      <td>-0.136953</td>\n",
       "      <td>-0.051486</td>\n",
       "      <td>-0.136953</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.575253</td>\n",
       "      <td>-0.575253</td>\n",
       "      <td>-0.144965</td>\n",
       "      <td>-0.575253</td>\n",
       "      <td>-0.144965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>-0.226133</td>\n",
       "      <td>-0.226133</td>\n",
       "      <td>-0.156852</td>\n",
       "      <td>-0.226133</td>\n",
       "      <td>-0.156852</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       pacf  pacf_ols   pacf_yw  pacf_ols_bymyself  pacf_yw_bymyself\n",
       "0  1.000000  1.000000  1.000000           1.000000          1.000000\n",
       "1  0.997827  0.997827  0.904116           0.997827          0.904116\n",
       "2 -0.700485 -0.700485 -0.122213          -0.700485         -0.122213\n",
       "3 -0.051486 -0.051486 -0.136953          -0.051486         -0.136953\n",
       "4 -0.575253 -0.575253 -0.144965          -0.575253         -0.144965\n",
       "5 -0.226133 -0.226133 -0.156852          -0.226133         -0.156852"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 使用statsmodels包的pacf函数\n",
    "from statsmodels.tsa.stattools import pacf,pacf_ols,pacf_yw\n",
    "from scipy.linalg import toeplitz\n",
    "\n",
    "# 自己实现用ols求解pacf的函数\n",
    "from numpy.dual import lstsq\n",
    "from statsmodels.tools import add_constant\n",
    "from statsmodels.tsa.tsatools import lagmat\n",
    "\n",
    "def cal_my_pacf_ols(x, nlags=5):\n",
    "    \"\"\"\n",
    "    自己实现pacf，原理使用的就是ols(最小二乘法)\n",
    "    :param x:\n",
    "    :param nlags:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    pacf = np.empty(nlags+1) * 0\n",
    "    pacf[0] = 1.0\n",
    "\n",
    "    xlags, x0 = lagmat(x, nlags, original=\"sep\")\n",
    "    xlags = add_constant(xlags)\n",
    "\n",
    "    for k in range(1, nlags+1):\n",
    "        params = lstsq(xlags[k:, :(k+1)], x0[k:], rcond=None)[0]\n",
    "        pacf[k] = params[-1]\n",
    "\n",
    "    return pacf\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def cal_my_yule_walker(x, nlags=5):\n",
    "    \"\"\"\n",
    "    自己实现yule_walker理论\n",
    "    :param x:\n",
    "    :param nlags:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    x = np.array(x, dtype=np.float64)\n",
    "    x -= x.mean()\n",
    "    n = x.shape[0]\n",
    "\n",
    "    r = np.zeros(shape=nlags+1, dtype=np.float64)\n",
    "    r[0] = (x ** 2).sum()/n\n",
    "\n",
    "    for k in range(1, nlags+1):\n",
    "        r[k] = (x[0:-k] * x[k:]).sum() / (n-k*1)\n",
    "\n",
    "\n",
    "    R = toeplitz(c=r[:-1])\n",
    "    result = np.linalg.solve(R, r[1:])\n",
    "    return result\n",
    "\n",
    "def cal_my_pacf_yw(x, nlags=5):\n",
    "    \"\"\"\n",
    "    自己通过yule_walker方法求出pacf的值\n",
    "    :param x:\n",
    "    :param nlags:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    pacf = np.empty(nlags+1) * 0\n",
    "    pacf[0] = 1.0\n",
    "    for k in range(1, nlags+1):\n",
    "        pacf[k] = cal_my_yule_walker(x,nlags=k)[-1]\n",
    "\n",
    "    return pacf\n",
    "\n",
    "cal_my_pacf_yw(x=data.values)\n",
    "\n",
    "# 结果比较\n",
    "\n",
    "pd.DataFrame({'pacf':pacf(data.value,nlags=5, method='ols'),\n",
    "              'pacf_ols':pacf_ols(data.value,nlags=5),\n",
    "              'pacf_yw':pacf_yw(data.value,nlags=5),\n",
    "              'pacf_ols_bymyself':cal_my_pacf_ols(x=data.values),\n",
    "              'pacf_yw_bymyself':cal_my_pacf_yw(x=data.values)})\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 偏自相关系数画图\n",
    "\n",
    "画图有两种方法：\n",
    "1. 直接使用`statsmodels`的`plot_pacf`;\n",
    "2. 使用`pacf`函数来计算数据，得到数据后再进行可视化，这种一般便于自己定制，也方便自己做更多的细节调整之类的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "(-1.0, 16.5)"
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x600 with 2 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAIHCAYAAACogAyRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAABV5UlEQVR4nO3deZwcdZn48c+T4ZAr4UgIgQQDHoirLAsIiwoEEcX1DpeC67murKjgsQKrIOD+FtEVQeMB6IIoyyFkQRclLmg4FKNBI7gKKiRrAklIICQQcsDk+f1R1dBpeqZnkunu6Z7P+/WqV6W+9a1vP10z0P3Mt+qpyEwkSZIkSX0b1e4AJEmSJGm4M3GSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkbZSI+PuIuCcinoyIR4dozIyIM4dirAavc2lEPN7s1+lEETEzImbWtI2PiGsi4uHyZ3Rye6KTpNbbpN0BSJI6V0S8CLgUuBH4HPBEm+N5MXAMcGlmzmtnLF3qS8BrgbOARcDs9oYjSa1j4iRJ2hhTKK5eOCkz/9zmWABeDHwGmAnMa2sk3elVwPWZ+e/tDkSSWs1L9SRJG2PHcv1oO4NQy+yIP2tJI5SJkyR1gYg4s7zn5EURcXVErCjvQ7kgIp5T0/c9EfGTiHgoItZExO8j4p/6GPd1EXFLRDxWjvmriDiu3DeP4pItgCWN7kuq3E8UEbtHxIyIWBkRD0bEGRERA3iPfxMRPyrjeDwibo6Iv63a/27ge+XmT8t4MiKmDGDsPmOKwryIuL7Occ+JiOURcWGD8TMipkXE0eX5XhURd0TES8v9H4iIP0fE6vLeoslVx55V3j82rs64F0XEo5WfcUTsV76PpeVrzI2I/6g5ZlREnBwR/1u+3uKIuDAitusn/ndHRAIBnFg5t/2eVEnqMiZOktRdrgaeA5wG/BD4CHBRTZ9/Av4P+Dfg48B84GsRcWJ1pzIRuQHYHjgHOBWYAxxRdjkZ+K+qMf8emN4gvh6K+6EWA58E7qRIvs7q76CI+CvgNuCvgc8DnwV2A2ZGxAFlt1uBL5f//rcynr8H/rAxMWVmAt8FXhcR29cc+0ZgdLm/kYOALwLfBs4E9gT+uzzvHwG+BnwBOBCoTna+Q3Fp/bHVg0XEZsBRwLWZuToidgR+DEymuN/sw8DlwN+yvgvL1/kZcBJwCXA8MCMiNu0j9lspziXA//DMuZWkkSMzXVxcXFw6fKH4Ip4U959Ut3+1bN+rqm2LOsffCNxXtT0GWAH8AnhOTd+o87pjBxDjpWXfL1ePBfw3sKZ6jLLfmVXb/1X22b2qbUIZ4y1VbUeVx04Z4HkbUEzAC8t+J9Qcfz0wt/qc9PE6CawGJle1/WPZvhDYpqr938r26r4/B35RM+Zbq98r8JZye79+4nhl2ee4mvbX1rZT3Cc2s877mNbu33cXFxeXdizOOElSd/lqzfZXyvXfVRoyc1Xl3xExJiLGArcAu0fEmHLX4cA2wOcyc3X1gJm5sZdoTasZaxqwGfDqep0jogd4DXBdZt5fdexC4D+BV0bE6GbGlJl/BGZRzMxU4toeeB1w+QDPyc25fqW/WeX62sx8rE777lVtlwEHRMTzqtqOp5gtvKXcfrRcv6GfmaOjgeXA/0TE2MpCMcv2OHDoAN6HJI1IJk6S1F3+VLN9H7CO4vItACLiFRFxU0SspPiyvYRilgOKmSaAyhf03w1xfOuA+2va/liuJ1PfOGBL4N46+/5A8Vk2qQUxXQa8IiKeW24fDWxKcSndQPylZnt5uZ7fR3v1PUdXUcyAHQ9Fwgu8gfWTtluAaymqCi6NiOvL+9k2rxrnBRQ/44cofu7Vy9Y8U+xDklTDxEmSutt6MyHljMXNwFjgY8DrKWaXvlR28XOhb1cCT/LMrNM7gNmZWS+hq6d3kO1PF8zIzGUUlw9WXvsoYHOq7q3KwlEU90hNA3ahuFfqzojYuuw2iiJpOryP5YwBvhdJGnF8jpMkdZcXUNxzU/F8ii/L88rtN1J84X5TZj49AxIRtZdo3VeuXwIM5fOZRlFcgvbHqrYXlut5z+pdWELxYN096ux7EcWMUWXWZkMuIxxQTJn5SETcABwfEZcDr6AokNEqlwHXR8TLKBKo32Tm/9Z2ysxfUNyb9qmyAuLlwNuAb1L8XF8N/Kz6kk1JUmP+ZVGSusuJNdsfLtc/KteV2Y2nZzPKy77eU3Pcj4HHgNPi2eXMG5YOb+BDNWN9iGIm5+Z6nTOzt4znzTVluscDxwG3Z+aKsnllud62STF9h+Ihu1+gOJdXDvJ1NsaPgKXAKcAh1FTyi4jt6vxs5pTryuV6V1NUETy9dvCI2CQith3CeCWpqzjjJEndZbeI+D5FlbwDKS4n+8/M/G25/8fAWuAH5bOHtgbeT3H51oTKIJm5IiI+SjFL8auI+E9gGUU58C2Bd21gfKuBIyLi2xRFEF5Hcbngv2Xmkn6O+zTFpWS3R8TXgKeAD1AkBJ+s6jeHIqE5pUwI1wA/ycyHhiimG4CHKe5v+lGDcYdUZj4ZEVdSJHW9wBU1Xd4FfDAi/otiZmkbip/tCorS9GTmLeXP/bSI2Jvi9+FJipnKoynKk1/T/HcjSZ3HGSdJ6i7HUiQLn6P48j8NeF9lZ3k/TqVk978DJ1A85+mC2oEy81vAmyi+eJ8OnAvswzOzVxuil+I5UDtRzNq8jOJ5Sc+aAamJ5X8pnoP0O4pnVH2G4llUh2bmrKp+i8r3tCPwLYrk4sVDFVNmrqUo1AADLwoxlC4r1zeXVQWr3QLMprgs78sUCeWfgFdl5tOXb2bmCRSl0HekKApyDvAqihmsnzU1eknqYLHxVWUlSe0WEWdSJBPjMnNpm8OpKyIuBY7KzK0b9R3OIuJLFMnoTpn5RItf+68pZtXemZntSNwkacRyxkmSpAEq7/d6B8Wzl1qaNJXeT/G8pelteG1JGtG8x0mSpAYiYkeKanRHATtQ59LGJr/+GykuOfxHYFpmrmxwiCRpiJk4SZLU2Ispyno/BHwkM+e0+PW/AoynKPLwmRa/tiQJ73GSJEmSpIa8x0mSJEmSGjBxkiRJkqQGRuQ9TuWT1XcGHmt3LJIkSZLabhvgweznPqYRmThRJE0L2h2EJEmSpGFjIvBAXztHauL0GMD8+fMZPXp0u2ORJEmS1CYrVqxg0qRJ0OBqtJGaOAEwevRoEydJkiRJDVkcQpIkSZIaMHGSJEmSpAaamjhFxMER8YOIeDAiMiLeMoBjpkTEryNiTUT8OSLeXafPiRExLyJWR8SsiNi/GfEPC729MHMmXHFFse7tbXdEkiRJ0ojT7BmnrYDfAicOpHNE7AbcAPwU2Bs4H/hmRLy2qs+xwHnAWcA+5fgzImLHoQx8WJg+HSZPhkMPheOOK9aTJxftkiRJklom+ilVPrQvFJHAWzPzun76nAu8PjNfUtV2JbBtZh5Rbs8CfpWZHyq3RwHzga9k5ucGGMtoYPny5cuHb3GI6dPhqKOg9ucTUayvuQamTm19XJIkSVIXWbFiBWPGjAEYk5kr+uo33O5xOhC4qaZtRtlORGwG7FvdJzPXldsHtijG5uvthZNOenbSBM+0nXyyl+1JkiRJLTLcEqedgMU1bYuB0RGxBTAW6Omjz059DRoRm0fE6MpC8WTg4eu222BBP8/nzYT584t+kiRJkppuuCVOzXIasLxq6ScrGQYWLhzafpIkSZI2ynBLnBYB42vaxgMrMnMVsBTo7aPPon7GPQcYU7VMHJJom2XChKHtJ0mSJGmjDLfE6Q7gsJq2w8t2MnMtcGd1n7I4xGGVPvVk5prMXFFZgMeGOvAhddBBMHHiM4UgakXApElFP0mSJElN1+znOG0dEXtHxN5l027l9q7l/nMi4rKqQ74B7B4Rn4+IF0XEB4FjgC9V9TkPeH9EvCsi9gS+TlH2/JJmvpeW6umBCy4o/l2bPFW2zz+/6CdJkiSp6Zo947Qf8JtygSLp+Q1wdrk9Adi10jkz5wKvp5hl+i3wceAfMnNGVZ+rgE+UY8yheN7TEZlZWzCis02dWpQc32WX9dsnTrQUuSRJktRiLXuO03DSEc9xqujtLarnffzj8MUvFpfnOdMkSZIkDYmBPsdpk9aFpA3S0wNTphQzT1OmtDsaSZIkaUQabsUhJEmSJGnYMXGSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkSZKkBkycJEmSJKkBEydJkiRJasDESZIkSZIa2KTdAagL9fbCbbfBwoUwYQIcdBD09LQ7KkmSJGmDmThpaE2fDiedBAsWPNM2cSJccAFMndq+uCRJkqSN4KV6GjrTp8NRR62fNAE88EDRPn16e+KSJEmSNpKJk4ZGb28x05T57H2VtpNPLvpJkiRJHcbESUPjttuePdNULRPmzy/6SZIkSR3GxElDY+HCoe0nSZIkDSMtSZwi4sSImBcRqyNiVkTs30/fmRGRdZYbqvpcWmf/ja14L+rDhAlD20+SJEkaRpqeOEXEscB5wFnAPsBvgRkRsWMfh0wFJlQtLwF6ge/V9Luxpt/bhzx4DdxBBxXV8yLq74+ASZOKfpIkSVKHacWM08eAizPzksz8PXAC8ATw3nqdM/ORzFxUWYDDy/61idOa6n6ZuayZb0IN9PQUJcfh2clTZfv8832ekyRJkjpSUxOniNgM2Be4qdKWmevK7QMHOMz7gCszc2VN+5SIeCgi7o2Ir0fEDv3EsXlEjK4swDaDeycakKlT4ZprYJdd1m+fOLFo9zlOkiRJ6lDNnnEaC/QAi2vaFwM7NTq4vBfqJcA3a3bdCLwTOAw4BTgE+FFE9DWdcRqwvGrpp/ybNsrUqTBvHvz0p7DPPsV67lyTJkmSJHW0TdodQAPvA+7OzF9WN2bmlVWbd0fEXcB9wBTg5jrjnENxn1XFNpg8NU9PD0yZUsw8TZnS7mgkSZKkjdbsGaelFIUdxte0jwcW9XdgRGwFvA34VqMXycz7y9d6fh/712TmisoCPDaA2CVJkiQJaHLilJlrgTspLqkDICJGldt3NDj8aGBz4LuNXiciJgI7AD4kSJIkSdKQa0VVvfOA90fEuyJiT+DrwFbAJQARcVlEnFPnuPcB12Xmw9WNEbF1RHwhIv42IiZHxGHA9cCfgRlNfSeSJEmSRqSm3+OUmVdFxDjgbIqCEHOAIzKzUjBiV2Bd9TERsQfwSuA1dYbsBfYC3gVsCzwI/Bg4PTPXNOEtSJIkSRrhWlIcIjOnAdP62DelTtu9QN0nqWbmKuC1QxmfJEmSJPWnFZfqSZIkSVJHM3GSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkSZKkBkycJEmSJKkBEydJkiRJasDESZIkSZIaMHGSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGTJwkSZIkqYGWJE4RcWJEzIuI1RExKyL276fvuyMia5bVNX0iIs6OiIURsSoiboqIFzT/nUiSJEkaiZqeOEXEscB5wFnAPsBvgRkRsWM/h60AJlQtz63Z/0ngI8AJwAHAynLM5wxt9JIkSZLUmhmnjwEXZ+Ylmfl7imTnCeC9/RyTmbmoallc2RERAZwM/GtmXp+ZdwHvBHYG3tKsNyFJkiRp5Gpq4hQRmwH7AjdV2jJzXbl9YD+Hbh0R/xcR8yPi+oj4q6p9uwE71Yy5HJjVYExJkiRJ2iDNnnEaC/QAi2vaF1MkP/XcSzEb9WbgHRQx/jwiJpb7K8cNeMyI2DwiRlcWYJtBvQtJkiRJI9qwq6qXmXdk5mWZOSczbwGmAkuAD2zEsKcBy6uWBRsfqSRJkqSRotmJ01KgFxhf0z4eWDSQATLzSeA3wPPLpspxgxnzHGBM1TKxj36SJEmS9CxNTZwycy1wJ3BYpS0iRpXbdwxkjIjoAV4KLCyb5lIkSNVjjqaorld3zMxck5krKgvw2ODfjSRJkqSRapMWvMZ5wLcjYjbwS4qKeFsBlwBExGXAA5l5Wrl9BvAL4M/AtsA/U5Qj/yYU5fYi4nzg0xHxJ4pE6rPAg8B1LXg/kiRJkkaYpidOmXlVRIwDzqYo3jAHOKKqxPiuwLqqQ7YDLi77LqOYsXp5Wcq84vMUyddFFMnV7eWY6z0oV5IkSZKGQitmnMjMacC0PvZNqdn+KPDRBuMlcEa5SJIkSVJTDbuqepIkSZI03Jg4SZIkSVIDJk6SJEmS1ICJkyRJkiQ1YOIkSZIkSQ2YOEmSJElSAyZOkiRJktSAiZMkSZIkNWDiJEmSJEkNbNLuACRJ0hDr7YXbboOFC2HCBDjoIOjpaXdUktTRTJwkSeom06fDSSfBggXPtE2cCBdcAFOnti8uSepwXqonSVK3mD4djjpq/aQJ4IEHivbp09sTlyR1ARMnSZK6QW9vMdOU+ex9lbaTTy76SZIGzcRJkqRucNttz55pqpYJ8+cX/SRJg2biJElSN1i4cGj7SZLWY+IkSVI3mDBhaPtJktZjVT1JUmtZKrs5DjqoqJ73wAP173OKKPYfdFDrY5OkLuCMkySpdaZPh8mT4dBD4bjjivXkyVZ7Gwo9PUXJcSiSpGqV7fPPN0mVpA1k4iRJag1LZTff1KlwzTWwyy7rt0+cWLT7HCdJ2mAmTpKk5rNUdutMnQrz5sFPfwr77FOs5841aZKkjdSSxCkiToyIeRGxOiJmRcT+/fR9f0TcFhHLyuWm2v4RcWlEZM1yY/PfiSRpg1gqu7V6emDKlGLmacoUL8+TpCHQ9MQpIo4FzgPOAvYBfgvMiIgd+zhkCnAFcChwIDAf+HFE1Fx3wI3AhKrl7UMevCRpaFgqW5LU4Vox4/Qx4OLMvCQzfw+cADwBvLde58w8PjO/lplzMvMe4B/KOA+r6bomMxdVLcua+SYkSRvBUtmSpA7X1MQpIjYD9gVuqrRl5rpy+8ABDrMlsCnwSE37lIh4KCLujYivR8QO/cSxeUSMrizANoN6I5KkjVMplV1b7a0iAiZNslS2JGnYavaM01igB1hc074Y2GmAY5wLPEhV8kVxmd47KWahTgEOAX4UEX1dxH0asLxq6edCe0nSkLNUtiSpww3rqnoRcSrwNuCtmbm60p6ZV2bm9zPz7sy8DngD8DKK+6PqOQcYU7VMbGbckqQ6LJUtSepgzU6clgK9wPia9vHAov4OjIhPAKcCr8nMu/rrm5n3l6/1/D72r8nMFZUFeGyA8UuShpKlsiVJHaqpiVNmrgXupKqwQ0RUCj3c0ddxEfFJ4HTgiMyc3eh1ImIisANgOSZJGu46sVR2by/MnAlXXFGsfd6UJI04m7TgNc4Dvh0Rs4FfAicDWwGXAETEZcADmXlauX0KcDZwHDAvIir3Qj2emY9HxNbAZ4BrKWatngd8HvgzMKMF70eSNJJMn148vLf6OVQTJxb3bDlTJkkjRtPvccrMq4BPUCRDc4C9KWaSKgUjdqV4DlPFPwGbAddQzCBVlk+U+3uBvYDvA38EvkUxq3VQZq5p4luRJI0006fDUUc9++G9DzxQtE+f3p64JEkt14oZJzJzGjCtj31TarYnNxhrFfDaoYpNkqS6enuLmabMZ+/LLKoBnnwyvPnNnXG5oSRpowzrqnqSpAHw/pvmuO22Z880VcuE+fOLfpKkrteSGSdJUpN4/03zLBxgvaGB9pMkdTRnnCSpU3n/TXNNmNC4z2D6qW/OmkrqAM44jXQzZzZv7Icfbu740kjW2wsf+EDf998AnHACjBkzfO+/Ge7/j+jthbFjYenSvvuMG/fMl/7haLifY4Bbb4WvfGX98zx2LHz4w3Dwwe2LS9LGmTKl3REMOWecJKkT3X13/1/oAZYsKfppw/T0FF/e+/OhDw3fxLQT3HorfOYzz/5dXrq0aL/11vbEJUl1mDhJUid6+OGh7af6Dj4YzjqrmAGpNm5c0e6MyIbr7S1mmvozbZqX7UkaNrxUT5I60Q47DG0/9e3gg+EVryhm7772NfjgB+GlL3WmaWMNZtZ0771bEpIk9ccZJ0nqRC996bNnQWqNG1f008br6Sm+vI8bV6xNmjaes6aSOoyJkyTV6u2FOXPg5puL9XC8VMj7b9TpnDVVN+mEzw1tNC/Vk6RqnVThq3L/TW2848YVSdNwi1eqVpk1bVS10FlTDXed9LmhjeKMkyRVdGKFr4MPhiuvhC99CV7wgmJ9xRV+WGv4c9ZU3aATPze0wUycJAk6u8KX99+oU1m1UJ2skz83tEG8VE+SwApfUrtYtVCdys+NEcfESZLACl9SO9XOmkqdwM+NEcdL9SQJrPAlSRocPzdGHBMnSQKfiyRJGhw/N0YcEydJAit8SZIGx8+NEcfESZIqrPAlSRoMPzdGFItDSFI1K3xJkgbDz40RoyUzThFxYkTMi4jVETErIvZv0P/oiLin7H93RPxdzf6IiLMjYmFErIqImyLiBc19F5JGDJ+LJEnt1dsLc+bAzTcX6+H+LCQ/N0aEps84RcSxwHnACcAs4GRgRkTskZkP1en/cuAK4DTgv4HjgOsiYp/M/F3Z7ZPAR4B3AXOBz5ZjvjgzVzf5LUmSJKlZbr21eLBs9TOSxo4t7ify0je1UStmnD4GXJyZl2Tm7ykSqCeA9/bR/yTgxsz8Qmb+ITNPB34NfAiK2SaK5OtfM/P6zLwLeCewM/CWpr4TSZIkNc+tt8JnPvPsB8suXVq033pre+LqRp02qzcMNHXGKSI2A/YFzqm0Zea6iLgJOLCPww6kmKGqNoNnkqLdgJ2Am6rGXB4Rs8pjrxyS4CVJktQ6vb3FTFN/pk0r7ifyUriN46zeBmn2jNNYoAdYXNO+mCL5qWenBv13qmob0JgRsXlEjK4swDYDiF2SJKk7dMLswt13P3umqdaSJUU/bThn9TZYZGbzBo/YGXgAeHlm3lHV/nngkMw8oM4xa4F3ZeYVVW0fBD6TmePLe6B+BuycmQur+lwNZGYeW2fMM4HP1LYvP+IIRm+66ca8xdZ54AHYZZehH/fhh4d+zIolS4qbJDuF8TZXp8ULnRez8TZfp8VsvM3VKfE+/DDMnQtr1z7TttlmsNtusMMO7Yur1pIl8Kc/Ne73ghcM3/M+3H8nMuHOO9f/Xai12Waw774QsXGvNZx+txpY8eSTjLnxRoAxmbmir37NLg6xFOgFxte0jwcW9XHMogb9F1W1LazpM6ePMc9h/cv/tgEWcNVVMHp0X7GPDDNntjuC4eNTn4L/9//aHcXAGa8ktVcn/H+tMrtQa+1auPfe4fWsoTlz4KMfbdzvgx8sKtdp8ObMgTvu6L/P2rVwzDEbf46nTNm441tpxQoYM6Zht6ZeqpeZa4E7gcMqbRExqtzu66d2R3X/0uFV/edSJE/VY44GDuhrzMxck5krKgvw2ODfjSRJUgcZ6D1Dw+WyvZe+9NkPkq01blzRTxtmoFcaNfOKpA7Wiqp65wHvj4h3RcSewNeBrYBLACLisog4p6r/BcAREfHxiHhReZndfsA0KK7FA84HPh0Rb4qIlwKXAQ8C17Xg/UiSJA1/nXbPUE9PUZygPx/6kIUhNsZAL5/roMvsWqnpiVNmXgV8Ajib4lK6vYEjMrNS3GFXYEJV/59TPLvpH4HfAkcBb6l6hhPA54GvABcBvwK2Lsf0GU6SJEnQmbMLBx9cXD5YO/M0btzwuqywUzmrt1Ga/gBcgMycRjljVGfflDpt3wO+1894CZxRLpIkSarVqbMLBx9clBy/++4iqdthh+KLvDNNG68yq1fvvrcKZ/X61JLESZIkSS1WmV3o73K94Tq70NNjAYhmqczq1T7Hady4ImlyVq9PJk6SJEndyNkF9cVZvQ3SiuIQkkayykMXlywZvg9dlKRu5T1D6ktlVu+ww4q1SVNDzjhJap5bb13/UoCPfrT48P7wh/2wlqRWcXZBGhImTpKao6+HLi5dWrT7l05Jah3vGZI2mpfqSRp6nfbQRUmSpAZMnCQNvU576KIkSVIDJk6Shl4nPnRRkiSpHyZOkoZepz50UZIkqQ8mTpKGXuWhi/0Zrg9dlKRGfMyCNCKZOEkaepWHLvbHhy5K6kS33gpve1vxeIU//alYv+1tRbukrmbiJKk5fOiipG5TecxCbfGbymMWTJ6kruZznCQ1jw9dlNQtBvqYhVe8wv/HSV3KxElSc/nQRUndYDCPWfD/eVJXMnEa6aZMaXcEw8cOO3TW+ei0eCWpky1cOLB+48f7/2apS3mPkyRJUiMTJgxtP0kdx8RJkiSpkYMOgokTIaL+/giYNKnoJ6krmThJkiQ10tMDF1xQ/Ls2eapsn3++hSGkLmbiJEmSNBBTp8I118Auu6zfPnFi0T51anviktQSFoeQJEkaqKlT4c1vhttuKwpGTJhQXJ7nTJPU9Zo24xQR20fE5RGxIiIejYhvRcTWDfp/JSLujYhVEfGXiPhyRIyp6Zd1lrc1631IkiStp6enqJz39rcXa5MmaURo5ozT5cAE4HBgU+AS4CLguD7671wunwB+DzwX+EbZdlRN3/cAN1ZtPzpUQUuSJElSraYkThGxJ3AE8LLMnF22fRj4YUR8IjMfrD0mM38HHFnVdF9EfAr4bkRskplPVe17NDMXNSN2SZIkSarVrEv1DqRIbmZXtd0ErAMOGMQ4Y4AVNUkTwFcjYmlE/DIi3hvRV23QQkRsHhGjKwuwzSBikCRJkjTCNetSvZ2Ah6obMvOpiHik3NdQRIwFTqe4vK/aGcBPgCeA1wBfA7YGvtzPcKcBnxlQ5JIkSZJUY1CJU0R8DjilQbc9Nzycp19nNHADxb1OZ1bvy8zPVm3+JiK2Av6Z/hOnc4Dzqra3ARZsbJySJEmSRobBzjh9Ebi0QZ/7gUXAjtWNEbEJsH25r08RsQ1F4YfHgLdm5pMNXm8WcHpEbJ6Za+p1KNuf3tfgyj5JkiRJWs+gEqfMXAIsadQvIu4Ato2IfTPzzrL5VRT3VM3q57jRwAyKJOdNmbl6AGHtDSzrK2mSJEmSpI3VlHucMvMPEXEjcHFEnEBRjnwacGWlol5E7ALcDLwzM39ZJk0/BrYE3gFUCjkALMnM3oh4IzAe+AWwmqLU+b8A/96M9yFJkiRJ0NznOB1PkSzdTFFN71rgI1X7NwX2oEiUAPbhmYp7f64ZazdgHvAkcCLwJSDKfh8DLh7y6CVJkiSp1LTEKTMfoe+H3ZKZ8yiSn8r2zOrtPo65kfUffCtJkiRJTdes5zhJkiRJUtcwcZIkSZKkBkycJEmSJKkBEydJkiRJasDESZIkSZIaMHGSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkjpNby/MnAkPPFCse3vbHZEkSVLXM3GSOsn06TB5Mhx6KPz618V68uSiXZIkSU1j4iR1iunT4aijYMGC9dsfeKBoN3mSJElqGhMnqRP09sJJJ0Hms/dV2k4+2cv2JEmSmsTESeoEt9327Jmmapkwf37RT5IkSUPOxEnqBAsXDm0/SZIkDYqJk9QJJkwY2n6SJEkaFBMnqRMcdBBMnAgR9fdHwKRJRT9JkiQNORMnqRP09MAFFxT/rk2eKtvnn1/0kyRJ0pAzcZI6xdSpcM01sMsu67dPnFi0T53anrgkSZJGgE3aHYCkQZg6Fd785qJ63sKFxT1NBx3kTJMkSVKTNW3GKSK2j4jLI2JFRDwaEd+KiK0bHDMzIrJm+UZNn10j4oaIeCIiHoqIL0SECaBGjp4emDIF3v72Ym3SJEmS1HTNTDguByYAhwObApcAFwHHNTjuYuCMqu0nKv+IiB7gBmAR8PJy/MuAJ4F/GarAJUmSJKlaUxKniNgTOAJ4WWbOLts+DPwwIj6RmQ/2c/gTmbmoj32vAV4MvDozFwNzIuJ04NyIODMz1w7h25AkSZIkoHmX6h0IPFpJmko3AeuAAxoce3xELI2I30XEORGxZc24d5dJU8UMYDTwV30NGBGbR8ToygJsM6h3o+7W2wszZ8IDDxTr3t52RyRJkqRhplmJ007AQ9UNmfkU8Ei5ry//CbwDOBQ4B/h74Ls14y6uOWZx1b6+nAYsr1oW9B++Rozp02HyZDj0UPj1r4v15MlFuyRJklQaVOIUEZ+rU7yhdnnRhgaTmRdl5ozMvDszLwfeCbw1Ip63oWOWzgHGVC0TN3I8dYPp0+Goo2BBTR79wANFu8mTJEmSSoO9x+mLwKUN+txPUbxhx+rGsvLd9uW+gZpVrp8P3Fceu39Nn/Hlus9xM3MNsKYqlkGEoK7U2wsnnQSZz96XWTxU9uSTi9LfVq2TJEka8QaVOGXmEmBJo34RcQewbUTsm5l3ls2vopjhmtX3kc+yd7leWK7vAD4VETtmZuVSwMOBFcDvBzGuRrrbbnv2TFO1TJg/v+g3ZUrLwpIkSdLw1JR7nDLzD8CNwMURsX9EvAKYBlxZqagXEbtExD0RsX+5/byIOD0i9o2IyRHxJopS47dm5l3l0D+mSJC+ExF/HRGvBf4V+Go5qyQNzMKFjfsMpp8kSZK6WtMegAscD9wD3Az8ELgd+Meq/ZsCewCVqnlrgVdTJEf3UFwWeC3wxsoBmdkLvAHopZh9+i5FclX93CepsQkThrafJEmSulpkvXs8ulxZknz58uXLGT16dLvDUTv09hbV8x54oP59ThEwcSLMnes9TpIkSV1sxYoVjBkzBmBMZq7oq18zZ5yk4aunBy64oPh3bbGQyvb555s0SZIkCTBx0kg2dSpccw3sssv67RMnFu1Tp7YnLkmSJA07XqrnpXrq7S2q5y1cWNzTdNBBzjRJkiSNEAO9VG+wz3GSuk9PjyXHJUmS1C8v1ZMkSZKkBkb0jNOKFX3OxEmSJEkaAQaaE4zUe5x2ARa0Ow5JkiRJw8bEzHygr50jNXEKYGfgsXbHMkDbUCR6E+mcmDuN57i5PL/N5zluPs9xc3l+m89z3Fye3+Zr5jneBngw+0mORuSleuUJ6TObHG7imecMPdZfpQ9tOM9xc3l+m89z3Hye4+by/Daf57i5PL/N1+Rz3HA8i0NIkiRJUgMmTpIkSZLUgIlTZ1gDnFWu1Rye4+by/Daf57j5PMfN5fltPs9xc3l+m6+t53hEFoeQJEmSpMFwxkmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkSZKkBkycJEmSJKkBEydJkiRJasDESZIkSZIaMHGSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSRoCI+PuIuCcinoyIR4dozIyIM4dirAavc2lEPN7s11FjG/ozj4jJ5bHvHvqoJKk1TJwkqctFxIuAS4H7gPcD/9jmeF4cEWdGxOR2xjFYEbFzGffeGzHG37Ui2ZQkDT0TJ0nqflMo/n9/UmZemplXtzmeFwOfASa3OY7B2pki7r03Yoy/K8eQJHUYEydJ6n47lutH2xmEJEmdzMRJklqkvMwrI+JFEXF1RKyIiIcj4oKIeE5N3/dExE8i4qGIWBMRv4+If+pj3NdFxC0R8Vg55q8i4rhy3zzgrLLrkkb3qFTuJ4qI3SNiRkSsjIgHI+KMiIgBvMe/iYgflXE8HhE3R8TfVu1/N/C9cvOnZTwZEVMGMHafMUVhXkRcX+e450TE8oi4sMH4h0fE7RHxaBn7vRHxb+W+KcCvyq6XVMX97nL/QRHxvYj4S/nzmh8RX4qILarGvxQ4sfx35fis2v+2iLiz6ud4d0ScVH3uymNeGRFfjoglZawXRsRmEbFtRFwWEcvK5fO1P7OI2CoivljGt6Z8j5+o02/zMv4lZTzfj4iJfZy3XSLiPyJicTnm/0bEe/s71+VxO0XEJRGxoDxuYURcHx12CaekkWOTdgcgSSPQ1cA84DTgb4GPANsB76zq80/A/wLfB54C3gh8LSJGZeZXK53KL+7/UfY9h2JW6W+AI4D/BE4ux31rOebjwF0N4usBbgR+AXyyHOssis+MM/o6KCL+CrgNWAF8HngS+AAwMyIOycxZwK3Al8v3/G/AH8rD//DsEQceU2ZmRHwX+GREbJ+Zj1Qd+0ZgNPDdBrH/N8W5OQNYAzwfeEVVfGcAZwMXle8T4Ofl+mhgS+DrwMPA/sCHgYnlPoALKS73Oxz4+5rXPxy4ArgZOKVs3rN8/Qtqwv0KsIjikr+/pbhn7VHg5cBfgH+huCTwn4HfAZeVrxEUv0+HAt8C5gCvBb4A7AJ8tOo1vgm8g+J36OfAq4Ab6py38RQ/kwSmAUuA1wHfiojRmXl+7TFVrgX+qnw/8yhmRg8Hdi23JWl4yUwXFxcXlxYswJkUXzCvr2n/atm+V1XbFnWOvxG4r2p7DEWS8gvgOTV9o87rjh1AjJeWfb9cPRZFUrGmeoyy35lV2/9V9tm9qm1CGeMtVW1HlcdOGeB5G1BMwAvLfifUHH89MLf6nNR5jZMbnSNgv7LPu+vsq/fzOhVYB+xa1Tat+Oh9Vt/zgeVATz+v/+7y9W+s+fn+vHydr1e19QDzgZlVbW8uj/9UzbjfK49/Xrn912W/r9b0u7zOz/ybwIPADjV9r6BI5rYotydXnztg23L7E63678/FxcVlYxcv1ZOk1vtqzfZXyvXfVRoyc1Xl3xExJiLGArcAu0fEmHLX4cA2wOcyc3X1gJmZbJxpNWNNAzYDXl2vc0T0AK8BrsvM+6uOXUgxa/HKiBjdzJgy84/ALOD4qri2p5gBubzBOXm0XL85Igb92Vjz89qq/Hn9nCLB+5sBDPEosBXFz7SRb9W8l1nl63yrKp5eYDawe1W/vwN6KWb8qn2xPP51Vf2o0+/86o1yButI4Afl5tjKAsygSOz36eM9rALWAlMiYrs++kjSsGLiJEmt96ea7fso/uI/udIQEa+IiJsiYiXFl+olFJe2QfGFFOB55fp3QxzfOuD+mrY/luvJ1DeO4lK1e+vs+wPF582kFsR0GfCKiHhuuX00sCnwnQbjXwX8jGIGZXFEXBkRxww0iYqIXaO4P+wRisshl1AkuvDMz6s/X6N4Pz8q7/n5j4g4oo++f6nZXl6u59dpr05Kngs8mJmP1fT7Q9X+ynodxe9ltdqf7TiKmaN/pHi/1cslZZ8dqSMz11Bckvg6ivN9a0R8MiJ2qtdfkoYDEydJar/1ZkIi4nkU97qMBT4GvJ5iJuJLZRf/3923KynurarMOr0DmJ2Z9RK6p5UzRgdTzF59B9iLIpn6n3I2rU/l/v+h+DmdC7yF4uf17rJLw59XZj5EUeb8TTxzH9KPIuLbdbr39jFMvfaGBT02QuV9fZfi/dZbftbXwVnc//RCinv9VgOfBf4QEQOZoZOklvPDV5Ja7wU128+n+P/xvHL7jcDmwJsy88LM/GFm3kRxeVO1yozAS4Y4vlGsf4kXFF9woe+b9pcATwB71Nn3IooZjMqMyIZcRjigmLIoCnEDcHw56/QKGs82VY5dl5k3Z+bHMvPFwKcoiiIc2iDul5axfDwzz83M68uf14P1Xqaf11+bmT/IzA9SzCZeCLwzIp4/kPgH4P+AnSNim5r2F1Xtr6xH8cyMZkXtz3YJ8BjFfVk39bE81F9AmXlfZn4xM19D8Xu8GfDxQb4vSWoJEydJar0Ta7Y/XK5/VK4rMwdPzxaU9zW9p+a4H1N8cT0tnl3OfGNnGj5UM9aHKGZybq7Xubyn5scU9whNrjp2PHAccHtmriibV5brbZsU03coHrL7BYpzeWWjgct7oWrNKdebl+u+4q738wrgJJ5tZbl/vTEiYofq7cxcxzPVDzdnaPyQomjEh2raP0qR0FV+/yrrj9T0O7kmxl6KynhHRsSzkveIGNdXIBGxZe3vLMUfAh5j6N6vJA0py5FLUuvtFhHfp6iOdiBl2efM/G25/8cUN87/IIpnD20NvB94iKJKHQCZuSIiPkpxX86vIuI/gWUUVdG2BN61gfGtBo4oLxObRXEfyuuBf8vMJf0c92mKy7Nuj4ivUZRR/wDFF+FPVvWbQ5FsnFImhGuAnzSYnRhMTDdQlAQ/GvhRo1mP0hkRcXB57P9R3JvzQWABcHvZ5z6K+81OiIjHKJKgWcA95b5/j4hdKKoIHsn69xdV3FmuvxwRM4DezLwS+GaZvP2kfM3nUiTUc2hcqn2gfgD8FPh/ZXL7W4qCHm8Gzs/M+wAyc05EXAF8sPz5/Bw4jGJmtNapFDNysyLiYuD3wPYURSFeXf67nhcCN0fE1eUxT1GUzB/PABJdSWqLdpf1c3FxcRkpC8+UBd+TogT0CuARiqp6teXE30jxxXYVRSntT1LMOCUwuU7fn1FcKrec4sv82+q87kDLkT9OcVncDIrkYFE5xqiavuuVpi7b/oYiIawkFj8BDqzzOv9AkWw8RYPS5IOJqeqYSon3tw/wZ/Mq4DrgAYpE7gGKaoAvqOn3JopnZj3J+uW196S4z+kxikvYLqK4T2q98uUUMz5fpkiC1/FMAcQjy/e2uHz9/wO+AexUdey7y/H26+P3amy981bTtjVwXvn+1lIUpPgENaXagedQPD9qaXnuv0/xTKp6P/MdKSoc/qUccyFwE/D+qj6Ta87XDuUxfyjHf5SirP7R7f7v1MXFxaWvJTI3tmKtJGkgIuJMioeWjsvMpW0Op66IuBQ4KjO3bncsGyMivgS8jyLxeKLd8UiSOp/3OEmSukp578w7gGtNmiRJQ8V7nCRJXSEidqS4r+YoikvBLmhvRJKkbmLiJEnqFi8GLqe4f+gjmTmnveFIkrqJ9zhJkiRJUgPe4yRJkiRJDZg4SZIkSVIDI/Iep/KJ7jtTPG9DkiRJ0si2DfBg9nMf04hMnCiSpgXtDkKSJEnSsDGR4gHhdY3UxOkxgPnz5zN69Oh2xyJJkiSpTVasWMGkSZOgwdVoIzVxAmD06NEmTpIkSZIasjiEJEmSJDVg4iRJkiRJDTQ1cYqIgyPiBxHxYERkRLxlAMdMiYhfR8SaiPhzRLy7Tp8TI2JeRKyOiFkRsX8z4pckSZIkaP49TlsBvwX+A5jeqHNE7AbcAHwDOB44DPhmRCzMzBlln2OB84ATgFnAycCMiNgjMx9qxptop7lLV3L17PksWLaKidttwTH7TWK3sVu1OyxJkiRpRIl+SpUP7QtFJPDWzLyunz7nAq/PzJdUtV0JbJuZR5Tbs4BfZeaHyu1RwHzgK5n5uQHGMhpYvnz58mFdHOLq2fM59dq7iAgy8+n1uUfuxdH7TWp3eJIkSVLHW7FiBWPGjAEYk5kr+uo33O5xOhC4qaZtRtlORGwG7FvdJzPXldsHtijGlpi7dCWnXnsX6xJ61+V661OuvYt5S1e2O0RJkiRpxBhuidNOwOKatsXA6IjYAhgL9PTRZ6e+Bo2IzSNidGWheDLwsHb17PlERN19EcFVs+e3OCJJkiRp5BpuiVOznAYsr1oWtDecxhYsW0Vfl1FmJguWrWpxRJIkSdLINdwSp0XA+Jq28cCKzFwFLAV6++izqJ9xzwHGVC0ThyTaJpq43Rb9zjhN3G6LFkckSZIkjVzDLXG6g6KSXrXDy3Yycy1wZ3WfsjjEYZU+9WTmmsxcUVmAx4Y68KF2zH6T+p1xOtbiEJIkSVLLNPs5TltHxN4RsXfZtFu5vWu5/5yIuKzqkG8Au0fE5yPiRRHxQeAY4EtVfc4D3h8R74qIPYGvU5Q9v6SZ76XVdhu7FeceuRejqiadeiIYFXDukXsx2ZLkkiRJUss0+zlO+wE/rdo+r1x/G3g3MAHYtbIzM+dGxOspEqWTKO5F+ofKM5zKPldFxDjgbIqCEHOAIzKztmBExzt6v0m8ZJfRvO6C2wF4zysn844DnmvSJEmSJLVYy57jNJx0ynOcAJ5Y+xQvPqPIG39/9mvZcrNm57qSJEnSyNGpz3GSJEmSpGHHxEmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkSZKkBkycJEmSJKkBEydJkiRJasDESZIkSZIaMHGSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGNml3AOo+c5eu5OrZ81mwbBUTt9uCY/abxG5jt2p3WJIkSdIGM3HSkLp69nxOvfYuIoLMJCK48Jb7OPfIvTh6v0ntDk+SJEnaIF6qpyEzd+lKTr32LtYl9K7L9danXHsX85aubHeIkiRJ0gYxcdKQuXr2fCKi7r6I4KrZ81sckSRJkjQ0TJw0ZBYsW0Vm1t2XmSxYtqrFEUmSJElDoyWJU0ScGBHzImJ1RMyKiP376TszIrLOckNVn0vr7L+xFe9FfZu43Rb9zjhN3G6LFkckSZIkDY2mJ04RcSxwHnAWsA/wW2BGROzYxyFTgQlVy0uAXuB7Nf1urOn39iEPXoNyzH6T+p1xOtbiEJIkSepQrZhx+hhwcWZekpm/B04AngDeW69zZj6SmYsqC3B42b82cVpT3S8zlzXzTaix3cZuxblH7sWoqkmnnghGBZx75F5MtiS5JEmSOlRTy5FHxGbAvsA5lbbMXBcRNwEHDnCY9wFXZmZtSbYpEfEQsAz4CfDpzHy4jzg2BzavatpmgK+tQTp6v0m8ZJfRvO6C2wF4zysn844DnmvSJEmSpI7W7BmnsUAPsLimfTGwU6ODy3uhXgJ8s2bXjcA7gcOAU4BDgB9FRE8fQ50GLK9aFgwwfm2A5+7wTJL0scNfaNIkSZKkjjfcH4D7PuDuzPxldWNmXlm1eXdE3AXcB0wBbq4zzjkU91lVbIPJkyRJkqQBavaM01KKwg7ja9rHA4v6OzAitgLeBnyr0Ytk5v3laz2/j/1rMnNFZQEeG0DskiRJkgQ0OXHKzLXAnRSX1AEQEaPK7TsaHH40xX1J3230OhExEdgBWLjBwUqSJElSH1pRVe884P0R8a6I2BP4OrAVcAlARFwWEefUOe59wHW1BR8iYuuI+EJE/G1ETI6Iw4DrgT8DM5r6TiRJkiSNSE2/xykzr4qIccDZFAUh5gBHZGalYMSuwLrqYyJiD+CVwGvqDNkL7AW8C9gWeBD4MXB6Zq5pwluQJEmSNMK1pDhEZk4DpvWxb0qdtnuBeHZvyMxVwGuHMj5JkiRJ6k8rLtWTJEmSpI5m4iRJkiRJDZg4SZIkSVIDJk6SJEmS1ICJkyRJkiQ1YOIkSZIkSQ2YOEmSJElSAyZOkiRJktSAiZMkSZIkNWDiJEmSJEkNmDhJkiRJUgMmTpIkSZLUgImTJEmSJDVg4iRJkiRJDZg4SZIkSVIDJk6SJEmS1ICJkyRJkiQ1YOIkSZIkSQ2YOEmSJElSAy1JnCLixIiYFxGrI2JWROzfT993R0TWLKtr+kREnB0RCyNiVUTcFBEvaP47kSRJkjQSNT1xiohjgfOAs4B9gN8CMyJix34OWwFMqFqeW7P/k8BHgBOAA4CV5ZjPGdroJUmSJKk1M04fAy7OzEsy8/cUyc4TwHv7OSYzc1HVsriyIyICOBn418y8PjPvAt4J7Ay8pVlvQpIkSdLI1dTEKSI2A/YFbqq0Zea6cvvAfg7dOiL+LyLmR8T1EfFXVft2A3aqGXM5MKvBmJIkSZK0QZo94zQW6AEW17Qvpkh+6rmXYjbqzcA7KGL8eURMLPdXjhvwmBGxeUSMrizANoN6F5IkSZJGtGFXVS8z78jMyzJzTmbeAkwFlgAf2IhhTwOWVy0LNj5SSZIkSSNFsxOnpUAvML6mfTywaCADZOaTwG+A55dNleMGM+Y5wJiqZWIf/SRJkiTpWZqaOGXmWuBO4LBKW0SMKrfvGMgYEdEDvBRYWDbNpUiQqsccTVFdr+6YmbkmM1dUFuCxwb8bSZIkSSPVJi14jfOAb0fEbOCXFBXxtgIuAYiIy4AHMvO0cvsM4BfAn4FtgX+mKEf+TSjK7UXE+cCnI+JPFInUZ4EHgeta8H4kSZIkjTBNT5wy86qIGAecTVG8YQ5wRFWJ8V2BdVWHbAdcXPZdRjFj9fKylHnF5ymSr4sokqvbyzHXe1CuJEmSJA2FVsw4kZnTgGl97JtSs/1R4KMNxkvgjHKRJEmSpKYadlX1JEmSJGm4MXGSJEmSpAZMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkSZKkBjZpdwCSJGlozV26kqtnz2fBslVM3G4LjtlvEruN3ardYUlSRzNxkiSpi1w9ez6nXnsXEUFmEhFceMt9nHvkXhy936R2hydJHctL9SRJ6hJzl67k1GvvYl1C77pcb33KtXcxb+nKdocoSR3LxEmSpC5x9ez5RETdfRHBVbPntzgiSeoeJk6SJHWJBctWkZl192UmC5atanFEktQ9TJwkSeoSE7fbot8Zp4nbbdHiiCSpe5g4SZLUJY7Zb1K/M07HWhxCkjaYVfUkSS1lqezm2W3sVpx75F6cUhaIAOiJIEnOPXIvJnueJWmDmThJklrGUtnNd/R+k3jJLqN53QW3A/CeV07mHQc816RJkjaSl+pJklrCUtmt89wdnkmSPnb4C02aJGkItCRxiogTI2JeRKyOiFkRsX8/fd8fEbdFxLJyuam2f0RcGhFZs9zY/HciSdpQlsqWJHWypidOEXEscB5wFrAP8FtgRkTs2MchU4ArgEOBA4H5wI8jYpeafjcCE6qWtw958JKkIWOpbElSJ2vFjNPHgIsz85LM/D1wAvAE8N56nTPz+Mz8WmbOycx7gH8o4zyspuuazFxUtSxr5puQJG0cS2VLkjpZUxOniNgM2Be4qdKWmevK7QMHOMyWwKbAIzXtUyLioYi4NyK+HhE79BPH5hExurIA2wzqjUiSNpqlsiVJnazZM05jgR5gcU37YmCnAY5xLvAgVckXxWV676SYhToFOAT4UUT09DHGacDyqmXBAF9bkjREKqWyR1VNOvVEMCqwVLYkadgb1uXII+JU4G3AlMxcXWnPzCurut0dEXcB91HcH3VznaHOobjPqmIbTJ4kqeUslS1J6lTNTpyWAr3A+Jr28cCi/g6MiE8ApwKvzsy7+uubmfdHxFLg+dRJnDJzDbCmauwBBS9JGnq1pbK33GxY/w1PkiSgyZfqZeZa4E6qCjtERKXQwx19HRcRnwROB47IzNmNXiciJgI7AAs3NmZJkiRJqtWKP/OdB3w7ImYDvwROBrYCLgGIiMuABzLztHL7FOBs4DhgXkRU7oV6PDMfj4itgc8A11LMWj0P+DzwZ2BGC96PJEmSpBGm6YlTZl4VEeMokqGdgDkUM0mVghG7AuuqDvknYDPgmpqhzgLOpLj0by/gXcC2FIUjfgycXl6SJ0nSkJq7dCVXz57PgmWrmLjdFhyz3yR2874sSRpRWnJheWZOA6b1sW9KzfbkBmOtAl47VLFJktSfq2fP59Rr7yIiyEwiggtvuY9zj9yLoy2hLkkjRisegCtJUkeau3Qlp157F+sSetfleutTrr2LeUtXtjtESVKLmDhJktSHq2fP77MSa0Rw1ez5LY5IktQu1oCVpA7n/TfNs2DZKjKz7r7MZMGyVS2OSJLULiZOI9wd9z085GOufrL36X/Puv8RnrNpz5C/hqTCzHsf4qLb7ieABAL4xi338YGDd+eQF+7Y5ujq66T/RzR66l/QnP+PbqxOOscAC5evYua9S1jy+BrGbb05U/YYx4QxW7Q7LEkb4cDn7dDuEIaciZMkdaiFy1dx0W33k1kkTfDM+sJb72eP8aPZacxz2hVeV5iyxzh+cNeDdfclcOgewzM57ST1kv8f3PXgsE7+JY1M3uMkSR1q5r1L+pwRCeCn9z7UynC60oQxW/CBg3en+janUQER8IGDdzcx3UjVyf+6ZL31hbfez6Llq9sdoiQ9zcRJkjrUksfXUP/um+Iv90se99F2Q+GQF+7IOW996dPbR7xkJ847em9nQ4aAyb+kTuKlepLUocZtvfnTlzfVinK/hsb40c/MLB2976Rhf89QpzD5l9RJnHGSpA41ZY9x/X7p9P4bDXeV5L8ek39Jw40zTpJUo1MqfFXuv7nw1uIeESjuv0m8/0adweIb6had8rmhjWPiJElVOq3C1yEv3JHJO2zFqdPvBor7bw7fcyeTJnUEk391g0773NCGM3GSpFKnlvf2/ht1MpN/dbJO/dzQhvEeJ0kqWeFLao/a5N8vmuoUfm6MLCZOklSywpckaTD83BhZTJwkqWSFL0nSYPi5MbKYOElSyfLekqTB8HNjZDFxkqRSpcJXVP35cFRAhBW+JEnP5ufGyGJVPUmqYoUvSdJg+Lkxcpg4SVINy3tLkgbDz42RoSWX6kXEiRExLyJWR8SsiNi/Qf+jI+Kesv/dEfF3NfsjIs6OiIURsSoiboqIFzT3XUiSJEkaqZqeOEXEscB5wFnAPsBvgRkRUfduuYh4OXAF8C3gb4DrgOsi4iVV3T4JfAQ4ATgAWFmO6ZyoJElSh1u4fBVX/PIvfPknf+KKX/6FhctXtTskqSWX6n0MuDgzLwGIiBOA1wPvBT5Xp/9JwI2Z+YVy+/SIOBz4EHBCRARwMvCvmXl9OeY7gcXAW4Arm/dWJEmS1Ewz732Ii267n6CoTBfAD+56kA8cvDuHvNAqdWqfpiZOEbEZsC9wTqUtM9dFxE3AgX0cdiDFDFW1GRRJEcBuwE7ATVVjLo+IWeWxA06cnlj7FJusfWqg3dviiar4nmhCrKuf7B3yMddUjbmmCeNLzdZpv8PG23ydFrPxqlMtWrGai267n0yeLvNdWV946/1M3mGr9e4nGi468Xd40YrV3PanJTz8+Fp22HozDnrBOHYawnPbjO+tzTLQWCOzr+rzGy8idgYeAF6emXdUtX8eOCQzD6hzzFrgXZl5RVXbB4HPZOb48lK+nwE7Z+bCqj5XA5mZx9YZc3Og+glk2wALJp18NaM233Kj36ckSZKkzrRuzRPMP/8YgDGZuaKvfiPlOU6nAcurlgXtDUeSJElSJ2n2PU5LgV5gfE37eGBRH8csatB/UVXbwpo+c/oY8xzWv/xvG2DBLz91GKNHj+4r9hFh1v2PtDuEtlvzZC8nXP5rAL5x/D5sPsxLiBqvJLVXp/1/7bY/LeGSn817+pK3UVFc/vbel+/GK18wtp2hPcv37pzPjb9bxLo6F0SNiuIZSUfvO6n1gXWRVp3jA3bffqPHaJUVK1Yw4fzG/ZqaOGXm2oi4EziMojoeETGq3J7Wx2F3lPvPr2o7vGwHmEuRPB1GmShFxGiK6npf7yOONcCaynaUj3fecrNN2HKzkf0oK58zsL7NN+3pqHNivJLUXsP9/2sLl6/ikp8/kzQBT39h/o+fz+Ulu4wZVg9qffWe4/nR7+r/bT2Bw/fcaVif706w7Ikn6etGnSz3D8U57qTv2E8NMNZWXKp3HvD+iHhXROxJkdxsBVSq7F0WEedU9b8AOCIiPh4RL4qIM4H9KBOtLG7KOh/4dES8KSJeClwGPEiZnEmSJAlm3ruE6GNfAD+996FWhtPQhDFb8IGDdyeimP2oXn/g4N2HVZLXqcZtvXm/vxPjtt68j71qeiqYmVdFxDjgbIpqeHOAIzJzcdllV2BdVf+fR8RxwL8C/wb8CXhLZv6uatjPUyRfFwHbAreXY65u7ruRJEnqHEseX9Pv7MKSx9f0sbd9DnnhjuwxfjQ/vfchljy+hnFbb86he+xo0jREpuwxjh/c9WDdfQkcuocl3/vSkjm0zJxGH5fmZeaUOm3fA77Xz3gJnFEukiRJqqMyu1AveRrOsws7jXkOb99/13aH0ZUqs3oX3rr+s7ISZ/Ua6ZyLDyVJkjQozi6oHmf1NoyJkyRJUpdydkF9cVZv8EycJEmSupizC9LQMHGS1FSLVjxTs+V7d87n1XuOZ8KYLdoYkSSNPM4uSBuvFeXIJY1QM+99iH/5r7uf3r7xd4v4+Pd+yy1/HF7lbyVJkhoxcZLUFAuXr+Ki2+4nq0o5rUvIhAtvvZ9Fy316gCRJ6hwmTpKaotMeuihJktQfEydJTdGJD12UJEnqi4mTpKaoPHSxnuH80EVJkqR6TJwkNcWUPcb1O+PkQxclSVInMXGS1BSVhy5GwKhgvbUPXZQkSZ3G5zhJahofuiipG/l8OmlkMnGS1FQ+dFFSN5l570NcdNv9T2/f+LtF/Oh3i/jAwbtzyAu9BFnqZl6qJ0mSNAA+n04a2ZxxGuEOfN4O7Q6h7Z5Y+9TT/z5g9+3ZcrPh/Z9Fp8UrSd3i3BvvYVQEvfns0jejIvjjQ4/x1n12aUNkklrBGSdJkqQBWLBsFVknaQLITBYsW9XiiCS1komTJEnSAEzcbgsi6j+hLiKYuJ0FIqRuZuIkSZI0AMfsN6nfGadj95vU4ogktZKJkyRJ0gDsNnYrzj1yL0YF9IyK9dbnHrkXk8du1e4QJTWRd5VLkiQN0NH7TeJlk7fnqtnzWbBsFRO324Jj95tk0iSNAE1LnCJie+ArwBuBdcC1wEmZ+Xg//c8CXgPsCiwBrgNOz8zlVf3qzZG/PTOvHNI3IEmSVMfksVtxyhEvancYklqsmTNOlwMTgMOBTYFLgIuA4/rov3O5fAL4PfBc4Btl21E1fd8D3Fi1/ehQBS1JkiRJtZqSOEXEnsARwMsyc3bZ9mHghxHxicx8sPaYzPwdcGRV030R8SnguxGxSWY+VbXv0cxc1IzYJUmSJKlWs4pDHEiR3MyuaruJ4pK9AwYxzhhgRU3SBPDViFgaEb+MiPdGX7VBSxGxeUSMrizANoOIQZIkSdII16xL9XYCHqpuyMynIuKRcl9DETEWOJ3i8r5qZwA/AZ6guB/qa8DWwJf7Ge404DMDilySJEmSagwqcYqIzwGnNOi254aH8/TrjAZuoLjX6czqfZn52arN30TEVsA/03/idA5wXtX2NsCCjY1TkiRJ0sgw2BmnLwKXNuhzP7AI2LG6MSI2AbYv9/UpIrahKPzwGPDWzHyywevNAk6PiM0zc029DmX70/saXNknSZIkSesZVOKUmUsoyoT3KyLuALaNiH0z886y+VUU91TN6ue40cAMiiTnTZm5egBh7Q0s6ytpkiRJkqSN1ZR7nDLzDxFxI3BxRJxAUY58GnBlpaJeROwC3Ay8MzN/WSZNPwa2BN4BVAo5ACzJzN6IeCMwHvgFsJqi1Pm/AP/ejPchSZIkSdDc5zgdT5Es3cwzD8D9SNX+TYE9KBIlgH14puLen2vG2g2YBzwJnAh8CYiy38eAi4c8ekmSJEkqNS1xysxH6Ptht2TmPIrkp7I9s3q7j2NuZP0H30qSJElS0zXrOU6SJEmS1DVMnCRJkiSpARMnSZIkSWrAxEmSJEmSGjBxkiRJkqQGTJwkSZIkqQETJ0mSJElqwMRJkiRJkhowcZIkSZKkBkycJEmSJKkBEyepw8x7eOXT/z7vf/7I3KUr++ktSZKkoWDiJHWQq2fP5w1fvv3p7Utun8dhX5zJ92bPb2NUkiRJ3c/ESeoQc5eu5NRr72JdPtPWm8m6hFOuvYt5zjxJkiQ1jYmT1CGunj2fiKi7LyK4ylknSZKkpjFxkjrEgmWryMy6+zKTBctWtTgiSZKkkcPESeoQE7fbot8Zp4nbbdHiiCRJkkYOEyepQxyz36R+Z5yO3W9SiyOSJEkaOUycpA6x29itOPfIvRgV0DMq1lufe+ReTB67VbtDlCRJ6lqbtDsASQN39H6TeNnk7blq9nwWLFvFxO224Nj9Jpk0SZIkNZmJk9RhJo/dilOOeFG7w5AkSRpRmnapXkRsHxGXR8SKiHg0Ir4VEVs3OGZmRGTN8o2aPrtGxA0R8UREPBQRX4gIE0BJkiRJTdPMhONyYAJwOLApcAlwEXBcg+MuBs6o2n6i8o+I6AFuABYBLy/Hvwx4EviXoQpckiRJkqo1JXGKiD2BI4CXZebssu3DwA8j4hOZ+WA/hz+RmYv62Pca4MXAqzNzMTAnIk4Hzo2IMzNz7RC+DUmSJEkCmnep3oHAo5WkqXQTsA44oMGxx0fE0oj4XUScExFb1ox7d5k0VcwARgN/1deAEbF5RIyuLMA2g3o3kiRJkka0ZiVOOwEPVTdk5lPAI+W+vvwn8A7gUOAc4O+B79aMu7jmmMVV+/pyGrC8alnQf/gaSeY9vPLpf5/3P39k7tKV/fSWJEnSSDSoxCkiPleneEPtssHlvjLzosyckZl3Z+blwDuBt0bE8zZ0zNI5wJiqZeJGjqcucfXs+bzhy7c/vX3J7fM47Isz+d7s+W2MSpIkScPNYO9x+iJwaYM+91MUb9ixurGsfLd9uW+gZpXr5wP3lcfuX9NnfLnuc9zMXAOsqYplECGoW81dupJTr72LdflMW28WG6dcexcvm7y9z0eSJEkSMMgZp8xckpn3NFjWAncA20bEvlWHv6p8vVl1B69v73K9sFzfAbw0IqqTssOBFcDvB/NepKtnz+8ziY4IrnLWSZIkSaWm3OOUmX8AbgQujoj9I+IVwDTgykpFvYjYJSLuiYj9y+3nRcTpEbFvREyOiDdRlBq/NTPvKof+MUWC9J2I+OuIeC3wr8BXy1klacAWLFtFZtbdl5ksWLaqxRFJkiRpuGraA3CB44F7gJuBHwK3A/9YtX9TYA+gUjVvLfBqiuToHorLAq8F3lg5IDN7gTcAvRSzT9+lSK6qn/skDcjE7bbod8Zp4nZbtDgiSZIkDVfR11/cu1lZknz58uXLGT16dLvDUZvMXbqSw744c717nCpGBfzk41O8x0mSJKnLrVixgjFjxgCMycwVffVr5oyTNKztNnYrzj1yL0YF9IyK9dbnHrmXSZMkSZKeNtiqelJXOXq/Sbxs8vZcNXs+C5atYuJ2W3DsfpNMmiRJkrQeL9XzUj1JkiRpxPJSPUmSJEkaIiZOkiRJktTAiL7HacWKPmfiJEmSJI0AA80JRuo9TrsAC9odhyRJkqRhY2JmPtDXzpGaOAWwM/BYu2MZoG0oEr2JdE7MncZz3Fye3+bzHDef57i5PL/N5zluLs9v8zXzHG8DPJj9JEcj8lK98oT0mU0ON0WeB8Bj/VX60IbzHDeX57f5PMfN5zluLs9v83mOm8vz23xNPscNx7M4hCRJkiQ1YOIkSZIkSQ2YOHWGNcBZ5VrN4TluLs9v83mOm89z3Fye3+bzHDeX57f52nqOR2RxCEmSJEkaDGecJEmSJKkBEydJkiRJasDESZIkSZIaMHGSJEmSpAZMnDpARJwYEfMiYnVEzIqI/dsdUzeIiNMi4lcR8VhEPBQR10XEHu2Oq5tFxKkRkRFxfrtj6SYRsUtEfDciHo6IVRFxd0Ts1+64ukFE9ETEZyNibnlu74uI06PqKYwanIg4OCJ+EBEPlv8/eEvN/oiIsyNiYXnOb4qIF7Qp3I7T3/mNiE0j4tzy/xEryz6XRcTObQy54zT6Ha7p+42yz8mti7DzDeQcR8SeEfH9iFhe/j7/KiJ2bWZcJk7DXEQcC5xHUXpxH+C3wIyI2LGtgXWHQ4CvAn8LHA5sCvw4IrZqa1RdKiJeBnwAuKvdsXSTiNgO+BnwJPA64MXAx4Fl7Yyri5wC/BPwIWDPcvuTwIfbGVSH24ris+zEPvZ/EvgIcAJwALCS4nPvOa0Jr+P1d363pPgu8dlyPRXYA/h+y6LrDo1+hwGIiLdSfMd4sBVBdZl+z3FEPA+4HbgHmALsRfF7vbqZQVmOfJiLiFnArzLzQ+X2KGA+8JXM/Fxbg+syETEOeAg4JDNvbXc83SQitgZ+DXwQ+DQwJzNPbmtQXSIiPge8IjMPancs3Sgi/htYnJnvq2q7FliVme9oX2TdISISeGtmXlduB8WXzC9m5r+XbWOAxcC7M/PKdsXaiWrPbx99Xgb8EnhuZv6lVbF1i77OcUTsAswCXgvcAJyfmee3PMAuUO8cR8SVwJOZ+fetjMUZp2EsIjYD9gVuqrRl5rpy+8B2xdXFxpTrR9oaRXf6KnBDZt7UsKcG603A7Ij4XnnJ6W8i4v3tDqqL/Bw4LCJeCBARfw28EvhRW6PqXrsBO7H+595yii+gfu41xxgggUfbHEfXKP/I/R3gC5n5v+2Op9uU5/f1wB8jYkb52Terv0smh4qJ0/A2Fuih+EtbtcUUHywaIuV/hOcDP8vM37U5nK4SEW+juCTktHbH0qV2p7iU7E8Uf9n8OvDliHhXW6PqHp8DrgTuiYgngd9Q/OX48vaG1bUqn21+7rVAefnjucAVmbmi3fF0kVOAp4AvtzuQLrUjsDVwKnAj8Brgv4DpEXFIM194k2YOLnWQrwIvofhLsoZIREwCLgAOz8ymXnc8go0CZmfmv5Tbv4mIl1DcH/Lt9oXVNY4BjgeOA/4X2Bs4PyIezEzPrzpWRGwKXA0ExR9fNAQiYl/gJGCf9H6YZqlM/FyfmV8q/z0nIl5O8dl3S7NfWMPTUqAXGF/TPh5Y1PpwulNETAPeAByamQvaHU+X2ZfiL0O/joinIuIpiqIcHym3e9obXldYCPy+pu0PQFMrC40gXwA+l5lXZubdmfkd4Es4g9oslc82P/eaqCppei7FH7acbRo6B1F87v2l6nPvucAXI2JeWyPrHkspZvRa/tln4jSMZeZa4E7gsEpbeUnZYcAd7YqrW5Qlb6cBbwVelZlz2x1TF7oZeCnFX+kry2zgcmDvzOxtV2Bd5GcUVbGqvRD4vzbE0o22BNbVtPXi52ezzKVIkKo/90ZTVNfzc28IVCVNLwBenZkPtzmkbvMdigpve1ctD1L8Eea17Qqqm5Tfj39FGz77vFRv+DsP+HZEzKaoenMyRYnGS9oZVJf4KsXlN28GHouIyvXzyzNzVfvC6h6Z+Riw3j1jEbESeNh7yYbMl4CfR8S/UHwZ2h/4x3LRxvsB8KmI+AvFpXp/A3wM+I+2RtXByiqbz69q2i0i9gYeycy/RPGct09HxJ8oEqnPUnzxvK7FoXak/s4vxQz1NRT3nb4B6Kn67Huk/EKqBhr9DgMP1/R/EliUmfe2LsrONoBz/AXgqoi4FfgpcATwRorS5M2Ly8svh7+I+BDwzxQ3xs4BPpKZs9oaVBcoy1vW857MvLSVsYwkETETy5EPqYh4A3AOxV+Q5wLnZebF7Y2qO0TENhRf3N9KcfnNg8AVwNl+ydwwETGF4otOrW9n5rvLkuRnUST/21I8q+WDmfnHVsXYyfo7v8CZFP+PqOfQzJzZlKC6TKPf4Tr952E58kEZyDmOiPdSXDY9EbgX+ExmXt/UuEycJEmSJKl/XqMtSZIkSQ2YOEmSJElSAyZOkiRJktSAiZMkSZIkNWDiJEmSJEkNmDhJkiRJUgMmTpIkSZLUgImTJEmSJDVg4iRJkiRJDZg4SZIkSVIDJk6SJEmS1ICJkyRJkiQ18P8BeHEams4p+6sAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sunspotarea = pd.read_csv(\"../datasets/sunspotarea.csv\",parse_dates=['date'])\n",
    "\n",
    "\n",
    "# 自己画一个\n",
    "# 自己计算\n",
    "pacf_value, pacf_interval = pacf(sunspotarea.value,nlags=15,alpha=0.05)\n",
    "\n",
    "xlabel = np.arange(start=0, stop=pacf_value.shape[0], dtype='float')\n",
    "\n",
    "fig, ax = plt.subplots(nrows=2, figsize=(10,6), sharex=True, dpi=100)\n",
    "ax[0].hlines(y=0, xmin=np.min(xlabel)-2, xmax=np.max(xlabel)+2, colors='red', linewidth=0.5)\n",
    "ax[0].scatter(x=xlabel, y=pacf_value, c='red')\n",
    "ax[0].vlines(x = xlabel, ymin=0, ymax=pacf_value, colors='red', linewidth=0.5)\n",
    "xlabel[1] -= 0.5\n",
    "xlabel[-1] += 0.5\n",
    "ax[0].fill_between(x=xlabel[1:], y1=pacf_interval[1:,0] - pacf_value[1:], y2=pacf_interval[1:, 1]-pacf_value[1:], alpha=0.25, linewidth=0, color='red')\n",
    "ax[0].set_title(\"pacf plot by myself\")\n",
    "\n",
    "\n",
    "# 使用别人写好的函数\n",
    "\n",
    "\n",
    "from statsmodels.graphics.tsaplots import plot_pacf\n",
    "plot_pacf(sunspotarea.value, ax=ax[1], lags=15)\n",
    "ax[1].set_title(\"pacf plot by statsmodels\")\n",
    "ax[1].set_xlim(-1, np.max(xlabel)+1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}